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ABSTRACT 



It is known that the static and dynamic behavior of 
fluids in porous media depends to a large measure on porous- 
media geometry*. In the past,, the ability to characterize 
this geometry has been restricted to such average properties 
as porosity and permeability. However, in recent years 
attempts have been made to achieve a more precise charac- 
terization based upon the fact that porous media is 
statistically composed. 

In this thesis techniques from statistical 'communi-' 
cation theory are adapted as possible methods for accomplish 
ing this classification process. Irt simplest form, a 
dichotamous function is defined by passing a line through a 
porous medium, the function having one value when the' line 
is in solid matrix, and another value when the line passes 
through pore space. The function is then analyzed using 
(l) Classical Fourier series harmonic analysis, and (2) 
determination of the autocovariance estimate and power 
spectrum. Comparisons are made between seyeral functions 
created from the same medium, and with functions created 
from other media. 

The results indicate that the autocovariance estimate * 
and the power spectrum, as characterizing functions, can 

ii i 



discriminate between different media. This success suggests 
many more paths of investigation, possibly leading to the 
complete characterization of porous media through statistical 
analysis . 
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CHAPTER I 



INTRODUCTION 

In August, 1965 > two American astronauts journeyed 
over two hundred miles into space in preparation for man's 
future journeys to the moon and thousands of miles beyond. 

As of August, 1965> man had journeyed only a few hundred 
feet beneath the surface of the earth upon which he lives. 

In this very interesting paradox, man is now capable of 
gathering information many miles into space, but has very 
little first-hand information about the sub-surface of the 
very earth from which he launches his space probes. If this 
seems incongrous, consider that while the benefits of space 
exploration are speculative at best, man derives the vast 

1 

majority of his resources, his energy sources, and even his 
food, from this earth. However, there are many explanations 
for .man's actions. Perhaps the foremost is the extreme' 
difficulty encountered in any attempt to explore extensively 
beneath the earth's surface. Costs would certainly be pro- 
hibitive, even if the engineering difficulties could be 
overcome. Another explanation might be that man has been 
able to theorize sufficiently about the earth's sub-surface 
areas from the limited information available so that he has 
not been willing or seen the necessity of risking great 
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expense, for what might be very little gain. A less prac- 
tical explanation, but in all probability a very real reason 
.for the lack of extensive sub-surface exploration, would be 
that there is less glamour to most people about traveling to 
the center of the earth than in a journey to Venus or Mars 
and thus, less interest in making the attempt. 

The above discussion is not meant as an implication that 
man is not interested in the earth's sub-surface. To the 
contrary, geologists, petroleum engineers, and others who 
derive their livelihood from the soil and beneath the soil, 
have performed extensive research in this area. From the 

exploration of # caves and mines, from core samples drawn from 

» • 

wells, and in many other ways, man has been able to learn 
much about the earth's constitution and composition. As a . 
result of this research, great improvements have been realized 
in the methods used to extract the natural resources from 
the earth, not to mention the improved conservation procedures 

t « 

which will allow man to benefit from the earth's resources 
for many years to come. 

As a part of this research, man has spent many years 
developing intelligence concerning the physical properties 
of porous media found in the strata of various kinds that 
lie beneath the earth's surface and the behavior of fluids 
within these porous media. Research efforts have been con- 
centrated on four properties which are considered fundamental 
and the basis for discussions regarding all other properties.^" 
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These four properties are porosity, a measure of the void 
space in the media; permeability, a measure of the media's 
ability to permit flow; the fluid saturation, and the 
electrical conductivity of the media. We will discuss these 
properties at greater length in Chapter II. 

The basis for this report is derived from a paper by 
2 

A. E. Scheidegger on "Statistical Geometry of Porous Media," 
published in 1961. As pointed out by Scheidegger, although 
the usual geometrical quantities, such as porosity and 
permeability, give generally good descriptions, there is a 
need for geometrical description of greater detail. This 
need arises because of the extreme complexity of the subject 
matter. For example, in some cases, media with a given 
porosity will hot allow as much flow as other media with less 
void space because of the lack of continuity in the pore 
spaces of the first sample. Media with identical permea- 
bilities to one fluid may have different permeabilities to 
another fluid. These are well-known phenomena in the area 
of petroleum engineering. As' a result, it is sometimes very 
difficult to classify media of different types without 
relying on extensive laboratory tests, which are expensive 
and time consuming. It is the aim of this report, to take 

advantage of .the inherent statistical nature of porous media 
1 » 

to characterize the media by a set of numbers or a mathe- 
matical function. To clarify the above statement, it is 
important to point out that porous media are consolidated in 
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a statistical nature. An individual particle has properties 
determined by its size, shape, density, electrical conduc- 
tivity, and various other characteristics. The total media 
will have properties determined by the properties of the 
individual particles, and by the means in which the individual 
particles’ have united to form the whole. This consolidation 
process will in turn depend upon the arrangement of the 
elementary particles, the relative number of different types 
of particles, and upon the relative number of large and- 
small particles. This dependency upon distribution functions 
is obviously statistical in nature. 

While much time has been devoted to research on the 

various physical properties mentioned above, very little has 

been accomplished in the latter area of investigating the 

3 4 

statistical nature of porous media. Scheideggar and Flood 
have provided impetus in this area, and defined some of the 
goals. It would be particularly desirable, for example, to 
be able to look at two pieces of porous media through a 
microscope and tell if they are identical in a statistical 
sense. It would be even more desirable to be able to make 
intelligent predictions about the other fundamental properties 
from the statistical description. 

The purpose of this report is to introduce and attempt 
techniques which utilize statistical communication theory to 
characterize porous media. In simplest form, the method will 
consist of defining a function created by passing a line over 
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or through a porous medium, and assigning a value of + l,when 
the line passes through pore space, and -1 when it passes 
through Solid matter. This function will then be analyzed 
in a number of ways. One attempt will be to fit the function 
with a single Fourier series and determine whether the 
Fourier coefficients effectively characterize the medium. 
Brief mention will be made of the double Fourier series 
whereby a surface area, rather than a single line, might be 
analyzed. 

Finally, the function will be looked at in terms of cer- 
tain statistical properties, derived from communication 
theory. The idea of random functions will be introduced and 
characterizing functions such as the autocovariance function 
and the power spectrum will be determined. Comparisons will 
be made within the same medium, and against different media. 
From the results, it is hoped to ascertain whether a given 
porous medium can be described and effectively classified by 
its power spectrum and autocovariance function, to the' 
exclusion of all others. 

t ■ 

If any of these three methods can be used to effectively 
classify porous media, the obvious step to be taken later 
will be the correlation of these descriptive quantities with 
the fundamental properties, such as porosity and permea- 
bility, which are now used to characterize porous media. 



CHAPTER II 



HISTORY 

Until quite recently, man has investigated and classi- 
fied porous media largely according to four fundamental 
properties from which it was felt all other important prop- 
erties could be derived. These four basic properties are 
permeability, porosity, fluid ’saturation , and electrical 
conductivity. Some of the earliest work was done by Henry 
Darcy, a Frenchman. In 1856 Darcy, investigating the flow 
of water through a sand filter, equated flow rate with the 
cross sectional area of the conduit, the length, and the 

t 4 

pressure head created by having the ends of the conduit at 

* I 

different levels. He found that a proportionality constant 
was required to balance the equality for various filters, and 
that the constant was a function of the sand packing. This 
constant became known as the permeability of the media. The 

I 

unit, of permeability is appropriately called the darcy, and 
is defined in the petroleum^ industry as: 

A porous medium has a permeability of one darcy 
when a single phase fluid of one centipoise viscosity 
that completely fills the voids of the medium will 
f-low through it under conditions of viscous flow at a 
rate of one cubic centimeter per second per square 
centimeter cross sectional area under a pressure or 
equivalent hydraulic gradient of one atmosphere per 
centimeter . 
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As can be seen from the .definition, Darcy's law may be 

extended to include other fluids as well as water. Ideally, 

the permeability is a function only of the porous medium 

and not the particular fluid flowing. However, investigators 

6 

such as Klinkenberg reported permeability variations dis- 
covered when using gases as the flowing fluid. The permea- 
bility was found to vary with the gas pressure, and the 
variation was attributed to a phenomenon known as gas slip- 
page. A discussion of this phenomenon is presented in most 
standard texts, e.g., Amyx, et al . ^ 

Prior to discussing the work of Carman and Kozeny in the 
prediction of the permeability of a porous media from its 
known pore size distribution, we should discuss early work in 

the area of calculating porosity. The early investigations of 
8 9 

Slichter and Fraser and Graton were concerned with the 
prediction of porosity of packings of uniform spheres. Since' 
porous media are, in most cases, quite obviously not composed 
of uniform spheres, but of a myriad of particle sizes and 
shapes consolidated with every conceivable angularity, the 
work of these men must be enlarged upon to deal correctly with 
real materials. Whereas their work was largely theoretical, 
and others continued to work with predictions of porosity based 
upon grain size distribution utilizing seive analysis, the 
more recent work in this area has been through empirical 
determinations of porosity, by laboratory measurement, using 
such devices as the Kobe porosimeter" 1 "® (mercury introduction), 
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or the Washburn-Bunting 11 porosimeter (air extracted by 
vacuum and measured ) . 

Studies of pore size distribution have been carried out 

12 

by several investigators including Carman, using gas 

13 

adsorption techniques; Ishkin and Kazaner, passing two 

14 

phase fluids through the media; Ritter and Erich, using 
small angle x-ray scattering patterns, and by many others 
with variety of techniques. Grain size distribution work 
has been accomplished mostly by simple seive analysis and 
there is apparently a need for more work in this area. In 
fact, Masson‘S indicates that after extensively researching 
this field, there is "a definite need for an inclusive 
mathematical function containing a few definite parameters 
that would define a given size distribution for most con- 
ditions approaching ideal grinding." 

One of the most interesting works done in the deter- 
mination of pore size distribution is attributed to Nuss and 
Whiting, , who devised a method of impregnating core samples 
with plastic under high pressure. The plastic was then 
polymerized at a high temperature and tfie rock material 
removed by acid leaching. The remaining plastic reproduces 
the original pore space. 

Many investigators have completed correlations between 

17 

the porosity and permeability of porous media. Carman, 
Kozeny, 1 ^ and Wyllie 1 ^ developed predictive equations from 
these correlations, whereby one property might be estimated 



9 



accurately from the other. In recent years in England, 

20 

Collis-George and Childs have refined the Kozeny Equation 

to improve results and allow greater generality. Francher, 

21 

Lewis and Barnes experimented with various porous matter 
to determine the relationship between grain size and fluid 
conductance. In their work they accounted for both laminar 

and turbulent flow. A brief discussion of the results is 

22 

found in the text by Amyx . 

Nearly all rock samples brought to the surface contain 
amounts of liquid entrained in their pore space. The deter- 
mination of the amounts and types of liquid is a very 
important characteristic description of the porous media. 

Comprehensive descriptions of the methods used to accomplish 

23 24 

this determination are covered by Hill, Horner, 

25 

Schilthius, and others. The two most common techniques 
for determination of fluid content are by retorting, and by 
solvent extraction and distillation. Fluid saturation may 
also be determined by using an indirect approach, that is, 
by measuring some other property and correlating the findings 
with past data to estimate the fluid entrained within the 
sample. Use of electric logs and capillary pressure measure- 
ments are two such indirect means of arriving at the fluid 
saturation . 

The fourth fundamental property of porous rock is the 
electrical conductivity. The rock, as it comes from the 
ground, is comprised of solid fragments, void -space and some 
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amount of entrained liquid. Generally, the solids and the 
oil and gas components filling the pore spaces are non- 
conductors of electricity. The water, if it contains dis- 
solved salts, then, is the only conductor. The ability or 
non-ability of a rock to conduct electricity is thus an 

indication of the interconnected or continuous water filled 

» « 

pores. Wyllie and Spangler created the first model of 
porous media taking advantage of this phenomenon to charac- 
terize the media according to its electrical resistivity, 
and hence, derived such other properties as tortuosity and 
porosity. Cornell and Katz 2 "' 7 and Wyllie and Gardner 2 ® 
added refinements in later models. 

From the f.our basic properties covered, many other 
properties can be derived. Researchers in many fields have 
developed various descriptive properties which characterize 
media more explicitly for their particular use. For example, 
the geologist has expanded the characterization to include 
mineral content, while the petroleum engineer, because of his 
interest in flow properties, has developed such concepts as 
capillary pressure, relative permeability, and wettability. 

A description of these properties is beyond the scope of this 
report, but the bibliography references the works in this 
area of Purcell, 2 ^ Leverett,®^ H. Brown, Rose ,® 2 Slobod,®® 

34 35 

Fatt, and Rapaport and Leas . ^ * 

Man has, then, taken advantage of many physical properties 
to descriptively classify porous media. However, in 1961, 
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A. E. Scheidegger and H. D. Fara, writing in the J ou rnal of 
Geophysical Research . point out that the usual geometrical 
properties, such as porosity and permeability, do not describe 
the media sufficiently for present requirements. Quoting 
their article: 

Geometrically, a porous medium is defined by 
giving the analytical equation of the surface which 
bounds its pore space. For any practical purpose, 
this is impossible to accomplish. Average geometrical 
quantities, such as porosity . . . , and specific 

surface . . . , are conceptually acceptable (and can 

be measured), but this is not true for the concept of 
"pore size." Since the pores form an extremely compli- 
cated system . . . , it is not easy to determine a 

pore size in a conceptually satisfactory fashion. 

Nevertheless, there is a great need for esta- 
blishing a better geometrical characterization of a 
porous medium than is afforded by the mere deter- 
mination of porosity and specific surface area. 

Scheidegger and Fara desired a means of characterizing porous 

media through analysis of the statistical make-up of the 

media. 

These investigators introduced four possible methods of 
statistical analysis: (l) by use of a correlation function, 

(2) a spectral analysis in terms of harmonic functions, (3) a 
spectral analysis in terms of other orthogonal functions, 
and (4) a spectral analysis of a specially constructed charac- 
teristic function of the media. 

In each of the four techniques, the writers developed 
statistical data for interpretation by drawing an arbitrary 
line through the porous medium. The line has a designated 
origin and data points are selected at equally spaced points 
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along the line. At each of these points, a value of + 1 is 
assigned if the line is, at that point, passing through void 
or pore space, and a value of -1 is assigned if the line is 
passing through solid. This series of 4 1 and -1 values, 
taken from the origin, define an orthogonal function f(s), 
which is characteristic of the line drawn through the function 
and of the media itself. If the medium is homogeneous and 
isotropic, and if the sample length is of sufficient size, 
theoretically, the functions or lines are independent of 
direction, origin, or position. The function is treated as 
a random function, as defined by Moyal,-^ and can be statis- 
tically analyzed in several ways. 

If the function is defined as above for void and solid 
spaces, the mean of the function may be derived: 

f - + 1(P) 4 (-1)(1-P) ' (2.1) 

where P ? porosity or pore spaces 

.’ . 1 = 2P ' 1 ( 2 . 2 ) 

In many cases, Scheidegger points out that it is incon- 
venient to work with a function having a mean other than zero, 
so he rewrites (2) as a new function with a mean that is 
equal to zero. The new function becomes: 

f’(s) - f ( s ) - (2P-1) (2.3) 

In the first of the four techniques- which the authors 
• * 

suggest, the function is characterized by its autocorrelation 



function: 
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R(A) = lim / + ® f(s) f ( s + A )ds lo M 

S >00 — . \ ) 

>: s s *■ 

Chapter III of this paper contains a discussion of the auto- 
correlation function. 

The second technique referred to in the Scheidegger 
•paper is to fit the function, f(s), with a Fourier series, 
and evaluate the Fourier coefficients, a^ and b^. A suitable 
evaluation suggested is to use the equation: 



c 



2 

k 





(2.5) 



2 

where c^ is no.t only a fair representation of the medium's 
geometry (statistically), but has the added advantage of being 
independent of point of origin or phase angle. 

The third and fourth techniques suggested are mere 
refinements of the two discussed briefly above. Each of the 
four techniques was offered b^ Scheidegger and Fara in an 
introductory manner, and there was very little in the way of 
tangible conclusions. However, these men opened an entirely 
new facet to the subject of characterizing porous media. 
Unfortunately, in correspondence between the author of the 
report in hand, and. Dr. Scheidegger, the latter indicated that 
he had done no further work in this area. Further works 'in 
the refinement of these and similar techniques are now left 
to this author and others. 

Perhaps brief mention should be included herein as to 
the history of statistical communication theory and time 



14 



series analysis, which Scheidegger has introduced into the 

o o 

world of porous matter. Granger^ credits the astronomers 
as being the first to analyze time series, although he men- 
tions discussions of meteorological and geophysical series 
in the middle of the nineteenth century. Even in these early 
discussions, it was noted that series of data taken over 
periods of time often follow patterns and contain certain 
trends, sometimes obvious and sometimes not so obvious. If 
the obvious trends such as those resulting from secular or 
seasonal variations could be subtracted from the series, the 
remaining series would contain random and independent data, 
and perhaps some of the not so obvious trends or patterned 
variations . 

Many of the works in separating these "hidden periodicities 

39 40 4l 

were instituted by La Grange, Buys-Ballot, Whittaker, 

42 43 

Stokes, and Schuster. J In all of these works, the random . 

independent series was assumed to be unimportant, and attempts' 

were made to mathematically generate similar functions as 

those creating the hidden trends. 

Schuster initiated the use of the periodogram, a method 

of pinpointing lesser periodic variations in a series. This 

approach was refined by such investigators as Whittaker and 

44 

Robinson who showed that apparent major trends may only be 
a combination of multiple lesser trends at their coincident 

peaks . • 

Zl ^ Zx * Zi/ 7 

It followed that Kolmogoroff, Weiner, and Cramer ' 

should refine these attempts at separating periodic tendencies 
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by using Fourier series analysis to segregate obvious sine 

or cosine trends. All of these earlier theoretical works 

led to practical applications, particularly in the field of 

48 

communications by Tukey and his assistants at the Bell 
Telephone Laboratories. Applications in other fields have 
followed, until today,' time series analysis and communication 
theory are major tools of the business and economic world 
as well as the world of science and engineering. 



CHAPTER III 



HARMONIC ANALYSIS OF RANDOM FUNCTIONS 

A simple means of analyzing a function is through a 
resolution into its harmonics by the use of a Fourier 
series. Where the function has one independent variable 
and one dependent variable, f(x) “ z, the one dimensional 
or single Fourier series is appropriate. A two dimen- 
sional function with two independent variables and the one 

dependent variable , • f ( x , y ) = z, requires the use of the 

# ♦ 

double Fourier series. A discussion of these mathematical 

tools is contained in Appendixes A and B. Where random 

processes or random functions are being considered, the 

classical harmonic analysis is not directly applicable. 

However, it has been shown that a "type" of Fourier analysis 
• « 

is still possible. This theory, which has .grown out of 
statistical theory of communication, provides the engineer 
with additional tools to assist in analyzing random 
phenomena. 

Two very important properties of a signal or function, 

i 

which for our purposes serve to characterize the function, 
are the power spectrum and the autocorrelation function. 
Prior to presenting the mathematical background for cal- 
culating these properties, it may be well to offer a lay- 
men's definition of each. 
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A. Power Spectrum 

This is a term used in basic communication theory which 

has been borrowed from the field of electricity. Various 

analogies have been drawn to describe the difficult concept 

of spectral analysis, one of the better ones being that 

49 • * 

discussed in Granger and Hatanaka. The author will 
attempt to draw an analogy of similar nature. If a person 
stands on the corner of 42nd Street and Broadway, during a 
very busy period, his ears will pick up many and varied 
sounds. Each of these sounds will appear to be different, 
although in many cases the "sound" which one hears is in 
reality a combination of many sounds coming from a range of 
signal creating devices. 











Figure, 3*1 The Human Ear Performing Simple 
Spectal Analysis 
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The different sounds are a result of many signal 
generators, emitting different frequencies. The human 
brain is equipped with a device, not unlike a filter used 
in electrical circuits, which isolates each signal and 
attempts to analyze it. Sometimes it is successful and 
reports,' "That is an automobile horn, and I must get out of 
the street," and sometimes the cacophony is so great that 
the individual signals can not be separated and interpreted. 
The brain is performing a type of spectral analysis when it 
selects, from a number of frequencies that the ear is hear- 
ing, individual signals and classifies them according to 
amplitude and frequency. Of course, in communication theory, 
and perhaps even in our analogy, there may be so many fre- 
quencies that it is only within our capability to identify 
and classify within ranges of frequencies. This is more 
often the case where a receiver is obtaining both intended 
or predetermined signals and "noise" which is added by some 
unknown and random generator. The mathematical process of 
decomposing a composite signal into its component frequencies 

and determining the proportion of the total signal at each 

* 

frequency has been called Power Spectrum Analysis, or it 
is often abbreviated as Spectral Analysis. 

Spectral analysis has been used more and more in 'recent 
years to analyze economic time series for the business 
world. Brown describes the use of this theory to separate 
the many factors which influence an economic series, such as 
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seasonal and secular trends. Here, as in electrical and 
communication work, we discover the presence of "noise," 
or factors which we cannot interpret as to origin, or 
predict as to future behavior. 

B. Autocorrelati on Function 

The autocorrelation function may be described as the. 
new function created when a function is correlated with 

I • 

itself. However, in order to provide a more meaningful and 
definitive .picture, we must look at the definition shown in 
Ekeblad . This text points out that the degree of auto- 
correlation in 'a time series indicates how well we may 
predict some as yet unknown future point in the series from 
the data available. A very simple analogy can be given. If 
the sales in a certain business fluctuate from month to 
month as shown in Figure 3-2. and we compare each month's 
sales with the preceding month's sales, we can see that fol- 
lowing every high sales month there is a lower sales month, 
and vice versa. 







— * » 1 1 — 
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Figure 3.2 Sales Function 
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It is very simple to predict that in the future for every 
high month, there will be a lower month to follow. However, 
if we increase our "lag" and compare only every second month, 
we find that it is harder to make such predictions. If we 
go further and compare only the third months, we lose even 
more confidence in our predictions. This is a very 
simplified explanation of Figure 3*3 which generally depicts 
the autocorrelation function. 




Lag in months 



Figure 3*3 Autocorrelation Function 

Figure 3-3 shows that' the autocorrelation function is 
maximum at a minimum .lag and generally decreases as the lag 
is increased. Relative maxima in the function appear which 
correspond to basic cycles existing in the data, but unless 
the function is periodic, i.e., repeats itself identically 
at some latter period, the value of the autocorrelation 
function will never approach the value determined .at zero 
lag. Likewise, as the lag increases, our confidence of 
prediction decreases. 
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While we will not make use of spectral analysis or the 
autocorrelation theory to separate and analyze signals, make 
predictions, or even use the correlation properties of the 
latter, we will attempt to take advantage of the properties 
of each to characterize or represent functions which are 
generated by our porous-media samples. For example, the 
autocorrelation function has a property of retaining all of 
the harmonics of a given function, containing no new ones, 
and discarding the phase angle. The usefulness of this 
property shall become apparent as we explain the application 
of these statistical methods to the work of characterizing 
porous media. 

With the background definitions of the power spectrum 
and the autocorrelation theory, we may now show the mathe- 
matical derivation of each. ' 

C . Periodic Functions 
<2 

Lee shows that a periodic function may be expressed 
as a Fourier series, as discussed in Appendix A. 



f (x) = + z a 

2 n =1 



cos + b sin 

L n 



nnx 

L 



(3-D 



where 




w \ unx J 

f ( x ) cos — : — dx 



(3.2) 
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o + L/2 
L J 

L -L/2 



f(x) 



sin 



nnx 

L 



dx 



' (3-3) 



L “ Fundamental period of the function, 
x “ value of the independent variable at equally 

spaced points along the function, such that x goes 
from -L/2 to L/2. 

n = number of harmonic terms to which the series is 
expanded . 

Equation (3*1) may be rewritten: 



f ( x ) - L F (n ) e jmi L (3-4) 

n = -°° 

in which: 

F(n ) =1 (a -jb ); forn =0, 1, 2, . . . (3-5) 

2 

j = the imaginary operator, equal to the square root 
of minus (negative) one. (Some texts use the 
notation "i," and a more definite description may 
be found in any text dealing with complex variables.) 

combining (3*2) and (3-3) 

F(n) - f / L//2 f (x ) e -^ nTT L dx . 

L -L/2 

*• 

Equation (3*6) is the Fourier Transform of the periodic 
function. Lee shows that this equation represents an 
"analysis" of the amplitudes and phase angles of each sinusoid 
into which function, f('x), is resolved. F(n) is also called 
the complex spectrum for this reason. 
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To separate the amplitude and phase characteristics of 



the complex spectrum, F(n), we may write it as: 

b 



F( 



v _ 1 / 2 2 

n) - — ♦b^ exp 



2 n n 



j tan 



-1 



n 



n 



where the amplitude is depicted by | F ( n ) j : 

| F ( n ) | ■ | X/V 

0 

and the term O n , indicates the phase angle: 



= tan 



-1 



b 

- _n 
a 

n 



(3-7) 



(3-8) 



(3-9) 



Combining the above, Equation (3*^) may be written: 



f(x) 



n =ro 

= ;s 

n = -°° 




cos 




(3.10) 



To analyze possible correlation between two functions, 



we use the function, 2 . 2 ’ w ^ ftre: 



X , 2 ) - r s . 

! , 2 L _ L / 2 1 



h -L/2 



f., (x ) f 2 (x + T) dx 



(3-11) 



We call the left-hand member of Eauation (3«H) the 
correlation function at a given lag ( 'T* ) . By multiplying one 
function by the other, at various values of Y , we can deter- 
mine the points of best correlation. 

If we adopt the special condition, where f-^(x) f 2 (x), 

we can write (3«H) as: 



24 






1,1 



{'Y ) 



i + L/2 

f I f, (x) f, (x + T ) dx 

L -L/2 1 1 



( 3 - 12 ) 



Lee shows that the autocorrelation function is also 
equat to: 



co 

r 

n = - c 



F(n) 



2 « Jntif T 



(3-13) 



When 'Y equals zero, Equation (3*13) reduces to: 



i * L / 2 ' 

t i i (T ) - i ; f ‘ 

1,1 L -L/2 1 



(x)dx 



CO 

I 

n *=-oo 



F 1 (n) 



( 3 - 14 ) 



or the mean square value of the function equals the sum, for 
all harmonics, of the square of the absolute value of the 
spectrum. The, power spectrum is then defined by: 



(n) 



F l(n) 



(3-15). 



It can be demonstrated that the power spectrum and the 
autocorrelation function are Fourier transforms of each other. 
As previously stated, all functions consisting of identical 
sinusoidal components have the same autocorrelation function, 
regardless of phase angle. Thus, the power spectrum, as a 
Fourier transform of the autocorrelation function, is like- 
wise independent of initial phase angle and having determined 
the latter, the former is uniquely determined. 

Lee also shows that these functions may be written in 
terms of Fourier coefficients: 
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If |F 1 (n)| 




then : 



4> (n) = | F( n )| 2 ° a n 




(3-16) 



and : 




cos ^ r (3.17) 



From Equations (3-16) and (3-17) , it is possible to 
calculate these desired properties of a function' with the use 



of a computer program. (See Appendix C.) 

D. Aperiodic Functions 

Whereas the tool of analysis for handling a periodic 

1 

function is the Fourier expansion, it is necessary to rely 

upon the Fourier integral to analyze the aperiodic function, 

The derivation of this integral and the steps necessary to 

evaluate the power spectrum and autocorrelation function for 

this type of function are not required within the scope of 

this paper. A discussion of the theory is contained in texts, 

53 54 

such as Tolstov ^ or Lee. 

E. Random Functions * 

We define a random function as follows: 

1. The generating process cannot be completely evaluated 
and a precise and unique determination of the value 
of the dependent variable for any value of the 
independent variable cannot be made. 
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2. As a conseauence of (1), the value of the dependent 
variable cannot be predicted at any given future 
point in time. 

In other words, a random function is one that is generated 
by some unknown source which dictates the manner in which the 
function behaves, but, due to its stachastic or probabilistic 
nature, cannot provide the actual value of the function at 
any given moment in time. An example of a random function 
might be the noise that a radio picks up in addition to the 
intended message or program. This noise may be the results 
of thermal emissions, static, or spillover from other fre- 
quencies, but an exact description of the noise is not pos- 
sible because of the complexity and random nature of the 
noise source. It might be noted that in most economic 
series, although containing obvious trends caused by seasonal 
or secular variations which may be isolated, certain strong 
random components exist. 

Because of the problems created by this "noise," it is 
of great import in communications work that analysis of such 
functions' be a possibility. In order to analyze random 
functions, we find it necessary to introduce certain re- 
strictions, and assumptions. First, we assume that the 
function has , an infinite period and th.at during this period, 
the function will trace a pattern, such that the value of 
the dependent variable will pass through, or infinitely 
close to, all possible values which are consistent with 
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the "ergotic hypothesis" as espoused by Maxwell, and 
discussed in Lee. Accepting this theorem, we may consider 
our random generating process as generating a group of 
finite length random functions rather than one infinitely 
long function having the combined characteristics of the 
entire grouping. The individual functions follow the nature 
of the infinite function on an average value basis and in a 
probabilistic manner. By the latter, it is meant that such 
properties as mean value crossings, or zero value crossings 
per unit length follow the probability distribution of the 
parent function. These finite functions will be referred to 
as an ensemble of random functions. They are all generated 
by the same source. This grouping of the ensemble is a very 
necessary step in the analysis of random functions, as 
realized and set forth by Gibbs. ^ 

In such an ensemble of functions, each of the finite 
functions will .be called member functions. It is important 
to realize that we assume that each of these member functions 
is generated by the same source, although we cannot mathe- 
matically formulate the operation of the source. This being 
so, we sense that properties which characterize any member 
of the ensemble, in turn characterize all others, as well as 
the infinite function from which the ensemble was created. 

One of these properties would be the autocorrelation function, 






rr ) 



lim 
L — ^ 00 



+ L 

L ■/ f (x)f n (x+T'Jdx 
2 -L 1 1 



1,1 



( 3 - 18 ) 
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which we assume exists for all values of 'Y . Another such 
property is the power spectrum. Wiener^ developed the 
theorem that the autocorrelation function of a random func- 
tion and the power density spectrum are related to each 
other as Fourier Cosine transforms. The fundamental assum- 
ption underlying this theorem is that- the autocorrelation 

function exists and is common for all member functions of 

57 

the ensemble. Blackman and Tukey state that, although 
Equation (3*18) is frequently called the autocorrelation 
function, a more appropriate name for the term defined by 
Equation (3*18) is the autocovariance function. They state 
that the autocorrelation function should be defined as the 
normalized autocovariance function obtained as the ratio 
of the latter function at any lag to the function at zero 
lag. This affords the autocorrelation function the property 
of having the value of unity at zero lag. However, it is 
convenient for our purposes, and entirely consistent with 
the definitions, to use the names interchangeably for 

t * 

characterizing porous media. 

In order that convenient equations might be developed 
which will allow us to calculate these properties with the 

c'O 

aid of a computer, we refer to Granger and Hatanaka? ’ They 

show that for. real processes the autocovariance '■'ir ’ may be 
• » 

written : 

Hr = 2 /q cos r u> f (uu)dtfj (3.19) 

where: tu ” < o 
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Formal inversion provides the power spectrum 



f(«0 

I- 




2 . cos 

j=z J 




( 3 . 20 ) 



For finite data, using numerical integration, an 
approximation to Equation (3*19) is: 



C T n-'T 



n-r 

2 

i-1 




( 3 . 21 ) 



where n = number of data points 

The estimate of Equation (3*20) becomes 
■ « 

A v / n-1 

f(u>) = ^ ( c 0 + 2 c j cos ^ 

6 0 

Hannan 0 shows that this is a poor estimate, but that 
it may be improved by properly weighting the values of the 
autocovariance terms, so that 



( 3 - 22 ) 
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(3.23) 



where "m" is arbitrarily selected by the user and de- 
notes the number of lags used. It should be less than 
one-third the value of "n." X ( ou ) represents the value of 



the weighting factors. 

There is considerable discussion concerning the choice . 
of the weighting factors, X(uj). However, the two most 
widely used, developed by J. W. Hamming and J. Von Hann , 
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are very similar and for this paper we will use the latter's 
system, now commonly known as the Tukey-Hanning estimate. 

It has the values 



2 [ ] 



COS 






m 



] 



( 3 -^ 4 ) 



The actual equations for estimating the covariance and 
the power spectrum may now be written 
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(3.25) 
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cos 4 
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( 3 . 26 ) 

i 



u j ■ L j-1 * -5 Lj * -25 L jn 

where * L 41 , 
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The L. terms are the raw estimate of the spectrum, and 

J 

the IK terms are the smoothed values. The smoothing 
operation of Eauation (3*27) is apparent from the form of 
the equation. 

To illustrate Equations (3*25) and (3*26), we might 
examine two very simple functions. Function A is an on-off 
function as shown in Figure 3*4. 
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Figure 3*^ Example Function A 

By examining the figure, it is seen that the function 
repeats every eight units and that the autocovariance 
function has the shape shown in Figure 3*5* 




Figure 3*5 Autocovariance Estimate of Function A 

Referring to Equation (3-26), it can be seen that the 
power spectrum will have a peak value when the frequency 
parameter "j” is such that the cosine term, cos , has a 

period coincident with the period of the autocovariance 
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function. As an example, taking m = 150 lags, the power 
spectrum reaches its peak value when the frequency parameter 
" j" is equal to 37* 5* The power spectrum appears as shown 
in Figure J.6. 



15 




Figure 3*6 Power Spectrum of Function A at 150 Lags 

These results are easily checked by observing that the 
original function passes through 3 7 n cycles in 150 units. 

A second function which can be used to demonstrate the 
nature of the analysis performed by use of the auto- 
correlation function and the power spectrum is Function B. 
Function B is a combination of two simple on-off functions 
which have been superimposed. The frequency of the first 
function is four cycles per eighty points. The frequency of 
the second function is ten cycles per eighty points. 

Function B, as shown in Figure 3 * 7 » has very little apparent 
regularity. 
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Figure 3*7 Example Function B 

The power spectrum for 50 lags shows peaks at 5 and 
12. 5j which were the individual peaks for each basic fre- 
quency used in the combined function. In other words, 
Equation ( 3 * 26 ) effectively selects the fundamental fre- 
quencies which constitute a function, even though they may 
not be apparent individually. Further, a measure of the 
relative contribution of each frequency is offered. 

The computer program "COVAR," discussed in Appendix C, 
has been written specifically to perform the mathematical 
operations set forth in Equations (3*25) and ( 3 * 26 ). The 
above figures were made from data and results of that pro- 
gram, run on Functions A and B. 



CHAPTER IV 



EXPERIMENTAL PROCEDURE 

A. General 

In ,this thesis an attempt is being made to classify 
or characterize porous media. It is our goal to differ- 
entiate between media which are not statistically identical 
through the use of vahious forms of statistical analysis. 

In order to apply the various tools outlined in Chapter III, 
we require a variety of samples which we know to be disimilar. 
In this chapter, we will discuss the selection and preparation 

of these samples and the collection of data therefrom. 

• * 

B. Preparation of Samples 

Data used .in this thesis were taken from two types of 
porous-media samples. The first sample was a photomicrograph 
of an actual porous medium. The other samples used were 
prepared by the author by imbedding multi-size spherical 
glass beads in plastic. 

The photomicrograph is a photograph of an actual porous 
medium taken through a microscope. In the sample used, the 
photographic image wa§ enlarged ninety-two times from an 
original size of approximately .12 centimeters by .16 
centimeters. This photograph was taken from the paper by 

Z Q 

Fara and Scheidegger. The original medium was sandstone 
with a porosity of 26.3 percent. a 
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The other, samples used were created by the author in the 
following manner: 

Materials : 

(1) Laboratory glass beads, three sizes: 

3mm . , 4mm., 5mm . 

(2) Natco Premium Quality Clear Polyester Casting 
Resin, NL-600. 

(3) Violet dye (water base). 

(4) Waxed paper cups. 

Procedure : 

.■ (l) Mix dye and resin to achieve a colored, 

slightly opaque plastic. 

(2) Place liquid plastic in waxed paper container. 

(3) Drop glass beads, in varying proportions as 

ft 

required, into the plastic. 

(4) Place 'in over for one day (30 degrees centi- 
grade). Allow to set. 

(5) Cut through sample in a random manner to 
realize an unbiased cross sectional surface 
of the mold. This surface area represents a 
model of a porous medium. Where the glass 
beads are apparent is taken as solid 'space, 
and where plastic is apparent is taken to be 



pore space. 
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For this paper, three samples were prepared. Sample X 
consisted of 3mm. beads only. Sample Y was a mixture of 
equal volumes of 4mm. and 5mm. beads. Sample Z was a mixture 
of equal volumes of 3mm., 4mm., and 5mm. beads. The cut 
samples were- approximately two inches in diameter. 

C . Preparation of Data 

Following ]the preparation of the samples, photographs 
of each were made by the Photo and Graphic Arts Department, 
University of Kansas (see Appendix D). The images of the 
created samples were expanded six times, while that of the 
Scheidegger photomicrograph was already of appropriate 
dimensions. An 80 x 80 grid, with ten units per inch was 
overlaid on each photograph. At each grid point, the data 
in point was recorded as being a 4 1 if a pore space (plastic), 
and a -1 if a, solid space (beads). These 80 x 80 data 
points were easily converted to input information on eighty 
column computer cards, for use with the programs written for 

% 

this thesis. (The negative values were read into the 
computer as fixed point zeros and changed internally to their 
correct values as negative ones.) 



CHAPTER V 



DISCUSSION OF RESULTS 

The mission of this project was to provide means of 
statistically analyzing porous-media, to show statistical 
variation in media known to be different, and to offer 
direction and inspiration to further work in this area. The 
author fully realizes that from the limited number of sam- 
ples dealt with in this thesis, there is insufficient evidence 
to prove that the far reaching goal of statistically charac- 
terizing all media can be accomplished. However, in this 
chapter, we will offer indications that 'the above outlined 
missions have been successfully completed and that success 
of the overall goal mentioned is entirely possible. It is 
hoped that from the results shown, others will be stimulated 
to continue research toward that end. 

Two basic computational tools were used to statisti- 
cally analyze the data taken from the porous-media samples. 
Although discussed in Chapter IV, it might be well to point 
out again that these data were taken at equally spaced points 
along lines randomly drawn across the cut-face of the samples. 
Where the line passed thrbugh pore space, a value of + 1 was 
assigned to the function, and where the line bisected solid, 
a -1 value was assigned. A series of these points taken in 

37 
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order along the line then formed a simple orthogonal function 
with the + 1 to -1 cycle arranged in a seemingly random pat- 
tern. Based upon the assumption that porous media are 
statistically constituted, it was hoped that several of these 
lines or functions taken from the same medium would be 
statistically similar while those taken from different media 
would vary by what is known as a "significant difference." 

The first method used to statistically analyze the 

functions was to fit each with a single Fourier series using 

"ID Fourier" computer program. (See Appendixes A and C.) 

Although each function was considered as being random in 

nature, it was analyzed over a finite length. The Fourier 

coefficients were determined and then combined according to 

2 2 2 

the Lee estimate of the power spectrum, c^ ” a n + b^ , so ' 
that the initial starting point or phase angle might be dis- 
regarded. These values of c n , as estimates of the power 
spectrum, are measurements of the contributions of the 
fundamentals frequencies constituting each function. This 
is essentially the technique used in the paper by Fara and 
Scheidegger . ^ 

Typical results of the single Fourier series approach 
are displayed in Table 5*1 Shown are the c^ values for four 
lines across the same media. As can be readily observed, 
there is very little consistency among the values, indi- 
cating that this method yields a poor estimate of the power 
spectrum. As a result, the simple Fourier approach was not 
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pursued to any great detail, although iater in this chapter 
we will compare the results of using this approach with those 
of our second technique on an identical sample. 

Because of the poor results achieved in the single 
Fourier approach, no analysis of the two dimension cross 
sectional faces of the media was attempted. However, a 
computer program "2D FOURIER" to determine double Fourier 
series coefficients for a -surface was developed, tested, and 
is fully operational. (See Appendixes B and C. ) 

The second statistical tool used to analyze the data 
was the determination of the communication theory properties 
of the functions or lines. These properties, the autoco- 
variance function and the power spectrum, are discussed in 
Chapter III. The properties were calculated using Fortran 
equivalents of Equations (3*25) and (3*26) in the computer 
program "COVAR." (See Appendix C.) As a part of that pro- 
gram, the power spectrum was plotted for each of several 
lines per media type, as well as the 90 and 95 percent 
confidence bands. These bands were created by multiplying 

the smoothed power spectrum by constants presented in Granger 
62 

and Hatanaka. The constant used is a function of the 
"equivalent degrees of freedom," an approximation chosen to 
allow for the fact that the series under consideration is 
not necessarily exactly normally distributed. These bands 
indicate 90 6r 95 percent confidence that the true value of 
the power spectrum lies within that range of values. In 
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other words, 90 or 95 percent, as the case may be, of all 
power spectra values from the same population should lie 
within that range. Generally, the approach was to compare 
data lines from the same media to see if the observed plots 
of the power spectra were similar. Then it was determined 
whether the values of each function within the sample, at 
each frequency, fell within the confidence bands of the 
other lines within that sample or type of media. A similar 
technique was used to compare the power spectra profile of 
one sample media with that of another media. 

For these comparisons, the effect of sample lengths was 
investigated. Samples of various lengths were created by 
combining the eighty point crossings developed in accordance 
with Chapter IV. The crossings were shuffled to reduce 
bias and eliminate peaks which would occur in the autoco- 
variance function at eighty point intervals. These peaks 
would be due to the similarity between adjacent lines drawn 
across a porous media. It was discovered after a series of 
trials that the longer the lines, the less bias there was 
to the idiosyncrasies of single crossings, i.e., a line 
might pass through continuous solid or pore space for an 
many, as twenty or thirty points without interruption or 
change. Longer lines tend to minimize the effect of these 
rare exceptions. As a result, lines of 3200 data points or 
40 crossings were utilized, this being the maximum number 
of points that could be handled on the computer at one time. 
It might be noted that lines of 1600 data points afforded 
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very similar results, indicating that the lines of 3200 points 
were of sufficient length. 

It was discovered that the number of lags selected, or 

more properly stated, the length of the maximum lag selected, 

played a very important part in the analysis. Again after 

many trials, it became apparent that when the maximum lag 

selected' was too short, the apparent frequencies were 

"packed" at the very low power spectrum values. In other 

words, with many frequencies constituting the function, 

♦ • 

each band or grouping of frequencies was too wide or contained- 
too many frequencies to effectively delineate the function. 

To long a maximum lag had the effect of causing random spikes 
in. the power spectrum which were characteristic of an in- 
dividual line 'but not of the media as a 1 whole. These "spikes" 

< « 

might be considered a result of sampling variability. Also, 
at long lags, values of the autocovariance function tended 
to distort the .power spectrum profile. In light of the 
above, a lag of 40 was selected as being optimum. It had 
the effect of smoothing the data and still offering good 
discrimination . 



spectra of different media. Appendix E contains the com- 
puter results from which these lines were plotted and the 
computer plots. The three profiles are obviously different, 
but it is important to discuss these differences in a more 
absolute statistical manner. As shown in the plots, the 




three plots of the power 
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95 percent confidence band of one of the lines is super- 
imposed upon the figure. By looking at each figure, it can 
be seen that generally the power spectrum of the other, or 
"outer" line lies within that confidence band. More 
specifically, the outer line is within the band at 24 points 
out of 40 in Figure 5 • 1 • This is the equivalent of 60 per 
cent of the time. In Figure 5*2 the result is 23 out of 40 
( 57 - 5 %) and in 5*3 it is again 23 out of 40 (57*5$) • 

ft 

While these results violate the 95 percent criteria, 
it is important to note that many of the points are very 
close to the limits and as only two examples were displayed 
per sample, variability of samples has a great effect. 
Ideally, many more lines per sample should be obtained, 
plotted and analyzed. Then the confidence bands around the 
mean of those lines should be superimposed upon a plot of 
the lines. Looking at all of the power spectra values that 
lie within these bands, the 95 per cent criteria would be 
approached. Using this method on four lines of shorter 
length from these same media, results very close to that 
criteria were realized and as more samples are considered, 
it is expected that the 95 percent criteria would be ful- 
filled. However, this was not proven in this work. By the 
time the tools to analyze media had been developed, and made 
available in the form of fully tested computer applications, 
and when the results discussed above became apparent, there 
was insufficient time to create or acquire additional samples 



for this thesis. 
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Cross comparing samples, there is very little correla- 
tion in comparison to the correlation within the same media. 
Between Figures 5*1 and 5 * 3 » which are very similar media 
(see Chapter IV), the base lines are within the other's 
confidence bands only 12 of 40 points, or 30 percent of 
the times. Between Figures 5.1 and 5«2, the result is only 
one out of 40 (2.5$), and between 5*2 and 5*3, it is one out 
of 40 (2.5$)« Visually comparing the different spectra, it 
can be seen that the 3mm* beads sample, which only has small 
beads in its make-up, has lower values than the other 
samples at low frequencies, and higher values than the others 
at higher frequencies. This is to be expected as there is 
a greater fluctuation in the 4 1 to -1 cycle in the original 
3 mm. bead data. There is less opportunity for long strings 
of 4 1 or -1 values due to this composition of the sample. 
Similar analysis can be made between the other samples. 

It is felt then, that this is a definite indication that 
the general shape of the power spectrum shown is indicative 
of the composition of the media from which the data is 
taken. Spectra from different media are not only visibly 
different, but also it is noted that the differences appear 
to be statistically significant. Again it is repeated that 
many more samples of various types need to be analyzed to 
prove this conclusively, but the technique is readily avail- 
able, and indications ■ are that the results would be 



favorable. 
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Of interest, Table 5*2 and Figure 5*^ display the com- 
puted autocovariance estimate and the power spectrum of the 
actual rock sample discussed in Chapter IV. 

A brief comparison was made between the two techniques 
used to analyze the functions, the single Fourier series and 
the power spectrum. Figure 5*5 shows this comparison. The 
function analyzed was the line drawn across the porous media 
shown in the Fara-Scheidegger paper. (See Chapter IV). As 
can readily be observed, there is very little consistency 
in the results. This is further indication that the single 
Fourier series approach is not particularly suited for this 
type of analysis. 
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FIGURE 5* j* 



POWER SPECTRUM ESTIMATION 
FOR 01SCRTTE T I HE SERIES 
USING TUKEY HANNING weighting FACTORS 
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TABLE 5 • 1 

Tabulated Results of Single Fourier Series 
Analysis of Three Size Beads Data 



Power Spectrum. 
Estimate 














1 


2 




3 


4 


C ( 0 ) * 


.2925 


.5637 




.4425 ' 


.2900 


C(l) 


.0005 


.0012 




.0048 


.0298 


C( 2 ) 


.0033 


.0002 




.0039 


.0029 


C ( 3 ) 


.0057 


.0000 




.0001 


.0029 


C(4) 


.0020 


.0002 




.0007 


.0006 


c( 5) 


.0078 


.0002 




.0017 


.0023 


C(6) 


.0002 


.0002 




.0017 


.0023 


c(7) 


.0028 


.0002 




.0023 


.0140 


C(8) 


.0058 


.0000 




.0041 


.0046 


C(9) 


.0022 


.0109 




.0018 


.0036 


C(10) 


.0027 


.0028 




.0017 


.0005 


C(ll) 


.0020 


.0011 




.0151 


.0071 


C ( 12 ) 


.0012 


.0016 




.0008 


.0018 


C(13) 


.0030 


.0000 




.0080 


.0050 


C(l4) 


.0005 


.0006 




.0040 


.0001 


c( 15) 


.0034 


.0003 




.0049 


.0018 


C(l6) 


.0049 


.0000 




.0012 


• .0024 


•C(17) 


.0048 


.0007 




.0017 


.0079 


C ( 18 ) 


.0049 


.0003 




.0007 


.0017 


C(19) 


.0202 


.0006 




.0087 


.0187 


C ( 20 ) 


.0003 


.0120 




.0048 


.0056 


C(30) 


.0024 


.0124 




* .0012 


.0015 


C (40 ) 


.0047 


.0139 




.0076 


.0039 


C(50) 


.0119 


. 00 66 




.0148 


.0010 


C(0) is an indication of the 
actual porosity of the sample 


porosity 
was 35 . 


of 

65 


the media, 
percent . 


The 
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SCHEIOEGGER DATA 3200 POINTS l AG * AO 



SMUOTHHD SPECTRUM WITH IN01CATE0 C0NF10ENCE LIMITS 









90 PCf. 




90 POT. 


9 5 PC T • 




95 PC T . 




COVARIANCE 


RAW 


IUWFR 




UPPER 


LOWER 




UPPER 


LAG 


SPFCTRUM 


SPECTRUM 


LIMIT 


AVERAGE 


LIMIT 


LIMIT 


AVERAGE 


LIMIT 


0 


0.0000 


0.653B 


0.64P2 


0.7171 


0.7716 


0 .6210 


0.7171 


0.7773 


l 


0.00*00 


0. 7003 


0.6429 


0.7112 


0.7652 


0.6159 


0.7112 


0. 7709 


? 


0.0000 


0.6303 


0.5622 


0.6219 


0 *6 69 l 


0.5305 


0.6219 


0.6741 


3 


o.ooco 


0.4465 


0.45R4 


0.5071 


0.5457 


0.4392 


0. 5071 


0.5497 


A 


0.0000 


0.5051 


0.4252 


0.4704 


0.5061 


0.4073 


0.4704 


0.5099 


5 


0.0000 


0. 42 A 7 


0.3890 


0.4303 


0. A630 


0.3726 


0.4303 


0.4664 


6 


0*0000 


0. 3666 


0.32 70 


0.3617 


0.3092 


0.3132 


0.3617 


0.3921 


7 


0.0000 


0. 2A09 


0.2643 


0.2923 


0. 3146 


0.2532 


0.2923 


0.3169 


0 


0.0000 


0.2249 


0.2055 


0. 2273 


0.2446 


0.1960 


0.2273 


0.2464 


9 


0.0000 


0.1705 


7). 1676 


0. 1853 


0. 1994 


0.1605 


0.1053 


0.2009 


10 


0.0000 


-0. 1755 


0.1443 


0.1596 


0.1710 


0.1302 


0. 1596 


0. 1730 


1 L 


0.0000 


0.1170 


* 0.1215 


0.1344 


0.1446 


0.1 164 


0.1344 


0. 1457 


12 


o.ooco 


0. 1201 


0.1061 


0. 1174 


0.1263 


0.1016 


0.1174 


0.1272 


13 


0.0000 


0.0962 


0.0P89 


0.0904 


0.1058 


0.0052 


0.0904 


0. 1066 


1 A 


0.0000 


0.0729 


0.0741 


0.0019 


0.0002 


0.0710 


0.0019 


0.0000 


15 


0.0000 


0.0557 


0.0691 


0.0764 


0.0823 


0.0662 


0.0764 


0.0829 


16 


o.ooco 


0.0614 


0.0635 


0.0703 


0.0756 


0.0600 


0.0703 


0.0762 


17 


0.0000 


0.0725 


0.0604 


0.0669 


0.0719 


0.0579 


0.0669 


0.0725 


IB 


0.0000 


0.061 1 


0.0565 


0.0625 


0.0673 


0.0542 


0.0625 


0.0670 


19 


0.0000 


0.0556 


0.0499 


0.0551 


0.059) 


0.0470 


0.0551 


0.0590 


20 


0.0000 


0 . 040 A 


0.04)3 


0.0479 


0.0516 


0.0415 


0.0479 


0.0520 


21 


0.0000 


0.0394 


0 .0302 


0.0423 


0.0455 


0.0366 


0.0423 


0.0450 


22 


0.0000 


0.0419 


0.0)00 


0.0420 


0.0452 


0.0364 


0.0420 


0.0455 


23 


0.0000 


0.0448 


0.0)70 


0.0410 


0.0441 


0.0355 


0.0410 


0.0444 


24 


0.0000 


0.032 3 


0.0323 


0.0350 


0.0305 


0.0310 


0.0350 


0.0300 


25 


0.0000 


0.0336 


0.0307 


0-03)9 


0.0)65 


0.0294 


0.0339 


0.0368 


26 


0.0000 


0.0362 


0.031) 


0.0347 


0.0373 


0.0300 


0.0347 


0.03.76 


27 


0.0000 


0.0327 


0.029) 


0.0325 


0.0349 


0.0281 


0.0325 


0.0352 


20 


0.0000 


0.0202 


C.0255 


0.0282 


0.0303 


0.0244 


0.0202 


0.0305 


29 


0.0000 


0.0235 


0.0227 


0.0251 


0.0270 


0.0217 


0.0251 


0.0272 


30 


0.0000 


0.0251 


0.0230 


0.0254 


0.0273 


0.0220 


0.0254 


0.0275 


31 


0.0000 


0.0279 


0.0243 


0.0269 


*0.0209 


0.0233 


0.0269 


0. 0291 


32 


0.0000 


0.0267 


0.0245 


0.0271 


0.0292 


0.0235 


0.0271 


0.0294 


33 


0.0000 


0.0273 


0.024) 


0.0268 


0.0209 


0.0232 


0.0260 


0.0291 


34 


0.0000 


0.0262 


0.0236 


0.0261 


0.0201 


0.0226 


0.0261 


0.0283 


35 


0.0000 


0.0248 


0.0215 


0.0238 


0.0256 


0.0206 


0.0230 


0.0250 


36 


0.0000 


0.0193 


0.0192 


0.0212 


0.0220 


0.0184 


0.0212 


0.0230 


37 


0.0000 


0.0213 


0.0205 


0.0226 


0.0244 


0.0196 


0.0226 


0.0246 


30 


0.0000 


0.0206 


0.0244 


0.0270 


0.0291 


0.0234 


0.0270 


0.0293 


39 


0.0000 


0.0296 


0*. 02 35 


0.0260 


0.0200 


0.0226 


0.0260 


0.0202 


40 


0.0000 


0.0163 


0.0200 


0.0230 


0.0247 


0.0199 


0.0230 


0.0249 



CHAPTER VI 



CONCLUSIONS AND RECOMMENDATIONS 

In this thesis, \an attempt has been made to provide 
effective methods of statistically analyzing porous media. 

It has been shown that the properties derived from communi 
cation theory, the autocovariance estimate and the power 
spectrum, are capable of differentiating between different 
media. Further, it would appear from the work with the 
simple representations of media which 'were created for this 
problem, that real porous media lends itself to this type 
of analysis. 

Perhaps the most significant conclusion that can be 

drawn from this thesis then, is that the techniques of 

statistical communication theory are adaptable for this form 

of analysis. Although the classical Fourier approach did 

63 

not provide consistent results, the Tukey Hanning power 
spectrum estimate offered good discrimination between the 
different media. 

There are several obvious steps required to further the 
proof and development of the theory presented herein. A 
large number of samples from each simple type of media should 
be analyzed to affirm the consistency of results. Then new 
samples should be fabricated with varying proportions of the 
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different size beads, to determine how fine a discrimination 
can be made. 

If these tools continue to offer good results with 
idealized media, the next step would be to test them on real 
rock samples. Such statistical methods as Discrimination 

Function Analysis or Factor Analysis might be applied to the 

» * 

data to improve the differentiation between "different" 
media . 

Finally, if these techniques prove successful upon real 

media, it would be logical to then attempt to relate the 

statistical properties, such as the power spectrum, to the 
• » 

static and dynamic fluid behavior of the media. 
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APPENDIXES 



APPENDIX A 



SINGLE FOURIER SERIES 

In many technical problems, it is necessary to deal with 
complicated periodic functions, that is functions F(x) = 

F(x + T) for all values of x, T being th’e period of the func- 
tion. Use of the Fourier series allows us to write these 
complex functions in terms of simple trignometric functions. 
Although Fourier series applications are not restricted to 
periodic functions, as explained in Gaskell, special care 
must be taken in applying the techniques, the discussion of 
which is beyond the scope of this paper. Most texts on the 
use of the Fourier series provide the necessary instruction 
in this area, e.g., Tolstov.^ 

Although the complete theory of Fourier series expansions 
are covered in detail in such reference material as Tolstov, 
Gaskell, ^ Hildebrand , ^7 and others, a resume of the simple 
derivation for a periodic orthogonal function is included 
herein to assist the reader in the basic understanding of the 
uses of the harmonic analysis that follow.* 

If we have a function, such that there is one independent 
variable and one dependent variable, taken between the limits 
-L < 0 < L, we could write the function: 

r + l 

z ** F(x) or; z m I xdx 

-L 



6l 
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We might illustrate the same function: 




,2 



-L 



0 



♦L 



x 

o 



x 



k 



. k 



The data set, (x^, z^) is such that x Q = -L, (x i + -^ - ) is 

constant and x^ “ + L, the value of "k" being even. 

We then assume that an equation using simple sine and 
cosine terms may be written, which will represent the func- 
tion, as follows: 

Z “ f(x) * + a^ cos ^ + a 2 cos + * * ' a n COS 



. -i . Tlx 4 , • 2TIX 

* b l Sln — + b 2 sin — 



b n sln 



nTTX 



or ; 

Z = f ( x ) = ^ + S (a n cos y + b n sin (1) 

n “1 

The degree to which this equation, so written, will converge 
upon the function, may be shown to be dependent upon two 
factors. The first factor is the number of harmonics (n) to 
which the series is expanded, and the second is the number of 
data points observed in relation to the basic length. The 
author has achieved an accurate fit to the degree of .97^ on 
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an eighty .point function of length 2 n , with data values 
ranging from -1 to *1 . The series was expanded to forty 
harmonics. By holding either of the stated factors constant, 
and decreasing ‘the other, a loss in accuracy of representation 
was realized. The reason for the "2" in the denominator of 
the first term of equation (l) will become apparent as we 
find the coefficients a n and b n of the series. By multiplying 
each side of the above equation by cos nx , and integrating 
between the limits of + L and -L, Gaskell shows that the 
equation may be rewritten: 



/ f(x) cos dx 

-L L 



a 



n 



;' L cos 2 dx 

-L ' L 



La 

n 



or: 



+L 

a = 1 / f(x) cos dx 

n L -L L 



( 2 ) 



With the a Q term as the first, term of the equation (1), 



utilizing the same steps, we discover: 

+L 



+ L 

/ f(x) cos 

-L 



( o )nx 
L 



+ L 



dx = / f(x) dx = / ~ La 0 

— L — L 



. + L , + L , . 

a = 7 / f(x) dx = 7 / f(x) cos \° dx 
° L -L L -L L 



( 3 ) 



a 



Gaskell further notes that is the average value of 

zrg 

f(x) over the entire interval . Using the trapezoidal rule to 
numerically integrate Equation (2) over its limits, and 



64 



2L 

utilizing equally spaced points (Ax = — ) where x = 1, 2, 3 > 
. . . k, we can approximate Eauation (2) as follows: 



n 






• i — !<■ n 

z 4 z. nTi . „ „ mix. 

o k cos 4 I z. cos i_ 

2 i “1 1 L 



2 

k 




nTi 



cos 



i=k-l 

4 £ 

i “1 



cos 



nnx^ 
~ L 



( 4 ) 



6 9 

It is shown by Whittaker and Robinson, that use of the 
trapezoidal rule produces coefficients which satisfy the 
least squares criterion of a "best" fit to the existing data, 
that is, 



i = k-l ? 

£ (z. - Px.) = Minimum (5) 

i =o 1 1 



where Px^ is the polynomial approximation to be observed z^. 

The same reasoning as above can be enlisted to determine 
the coefficient b n , except that the sine terms are zero at 
the end points, (sine (0) = 0) and b n becomes: 



n 



, 1-k-l 

= I S 

K i-1 



. mix. 
z^ sin l 



( 6 ) 



Now, assuming we have discrete data z for each value 
of x^ taken at equally spaced points between the limits, we 
can now compute the values of the coefficients, and having 
so done, substitute back into Equation (1), which becomes a 
one-equation representation of the function. The coefficients, 
a n and b , serve to characterize, in one way, the function. 
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However, this means of characterization has the distinct 

disadvantage of being dependent upon the phase angle of the 

70 

sinusoids into which the function is resolved. This 

disadvantage may be overcome in several ways. Some of the 
more complex characterization techniques, which we shall 
discuss and compare in latter sections, are the uSe of the 
autocorrelation function, and the power spectrum. These 
functions are non-dependent of phase angle, as we will 
later show. A less complex technique, which still dismisses 
the dependency upon the phase angle is discussed in a paper 
on wind generated sea waves, distributed by New York 
University .7^ An estimate, of the power spectrum is made 
in the following manner, using the Fourier coefficients: 

2 2 2 

Estimated Power Spectrum c^ = a n + ^*n (7) 

Although this estimate is a characterization of the 

function, and does dismiss the dependency upon the phase 

2 

angle, the same reference notes that the quantity c n is a 
very unstable estimate of the energy at any frequency. : Only 
through a weighting process of the c n values, would an 
accurate spectrum, such as that determined by the Tukey 
method"^ be recovered. Lee^ uses a similar estimate of the 
power spectrum, with the exception that he divides the above 
estimate by a factor of four. 



n 



n 



n 



4 



( 8 ) 



In that the basis for our discussion dealing with auto- 
correlation and power spectrum analysis evolves from this 
same reference, and since the division does not affect the 
term for our purposes, we have used the "Lee" estimate in 
our work. This usage will allow a means of comparing this 
simple estimate of the power spectrum with the values 
computed by more complex techniques. 
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DOUBLE FOURIER SERIES 



As in the case of a function with one independent 
variable discussed in an earlier section, it is also possi- 
ble to represent with a Fourier series, a dependent 
variable as the function of two independent variables. 

z = f(x,y) (1) 

A rigorous discussion regarding orthogonal systems in two 
variables' may be found in Tolstov.?^ 

Thus, a parameter such as z may be considered to be 

« 

oscillatory in two mutually perpendicular directions, such 
as east-west and north-s’outh, allowing Equation (1) to ■ 
represent the areal variation in z . If we consider a 
function z to have a fundamental period of 2L along the x 
coordinate, and 2H along the y axis, then following the same 
reasoning discussed earlier in the single dimension case, we 
can write Double Fourier series function: 



m=M n = N 

z “ p = r f x 

m=o n=o m ’ n 



a cos TimX cos Tiny 

m ’ n IT H 



b sin Timx cos Tiny 4 c cos nmx sin Tiny h 
m ’ n L h m,n L H 



d 



m , n 



sin Timx sin nny 
L . H 



( 2 ) 



6? 



where 






1/4 


o 

II 

c 

11 

B 


m,n' 1/2 


m = o,n>o or m>o,n = o 


“ 1 ' 


m>o , n>o 


m,n 


• 


» N 


number of harmonics < 



(Should be equal to or less than one 
half the number of data points con- 
sidered in the x and y directions 
respectively.) 75 

The reasoning behind the quarter and half terms follows also 
from the one dimension analysis, in that the end terms con- 
sidered in a numerical integration are averaged, according 
to the theopy of the trapezoidal rule. 7^ in the two dimen- 
sion case, the averaging process must be carried in each 
direction, thus the quarter terms. Again, the problem is to 
evaluate the coefficients. Integrating Eauation (2), the 
coefficients are determinable from the following relations: 

, + H 4 L 

a - f „ ; / 

m,n LH ^ 



f(x,y) cos 



=£* 003 aydx 



(3) 



m,n 



m ,n 



m,n 



h’ 



h' 



h> 



4 H + L 

; 

-H -L 

4 H + L 

I 

-H -L 

♦H + L 
J 

-H _L 



f(x,y) sin — P- cos dydx 



f ( x , y ) cos 



f (x,y) sin 



™ sin HSZ dydx 



Timx 



sin 



Tiny 

H 



dydx 



(4) 



(5) 



-H 



( 6 ) 
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It is possible to numerically integrate Eauations 3-6 
utilizing the trapezoidal rule. Again it must be remembered 
that the integration must be carried out in two dimensions, 
and special care must be taken with regard to the corner 
points. In order to clarify the following derivation, 

Figure 1 displays the nomenclature to be used. 



y 




j=o,l, 

i*o,l, 



2 , 

2 , 



. X 
. k 



Figure 1 

According to the proper use of the trapezoidal rule, 
it can be shown that- for the one dimension case (one 
independent variable) 

+ L 

/ z(x) dx = Ax 

-L 

where Ax = ^ 

Considering the second dimension of depth in our picture, 
when we are concerned with but one direction, we are able to 
write Equation (7) 



z ♦ z i=k - 1 
fo fk * L z. 

2 i-1 



(7) 
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. + L 

f z ( x ) dx 

-L 





i-k-1 

♦ Z 

i-1 



z 



i, J 



( 8 ) 



However, in a specific case, where it is required that 
integration be applied in two dimensions, as will be done for 
Equation (3)> the following steps are taken: 



a 



m ,n 



1 1 
L * H 



LH -H 



+H 

'-H 


+ L 

/ 

-L 


z(x,y) 


dxdy 


Ax 


z 

' O, ,1 


+ 

N 


i=k-l 
♦ Z 
i-1 




2 




dy 



_ Ax Ay 
~ LH 





+ 



2 



j=^-l 

Z 


z . 4 z, 

O, .1 k , .i 


(B) 

+ 


1=k_1 Z + Z n 

z i.o z iJ 


n 


2 




i=l ,.2._ 


- 









j = £-l 

Z 

n 

where: (A) 



(B) 



(C) 



(D) 



i=k-l 

Z 

i-1 



'i, J 




(9) 



term that averages the corners in both 
directions. 

term that averages opposite vertical edges 
( less corners ) . 

Term that averages opposite horizontal 

edges (less corners). 

term that sums all interior points. 
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The number of points in either direction should be an 
odd number, requiring "k" and " " to bS even numbers. Carry- 
ing the procedure to completion and reinserting the trigno- 
metric terms, which were omitted, to simplify the explanation, 
we can write the coefficient as follows: 



a - 1 . §£ . m. 

m,n k 1 





z^ + z, z 

0,0 k ,o + o, 


1 

#> 

N 

+ 




2 


2 cos mTT cos nn 




2 





r*-i 

z 

j=l 



Z z o , ,j 4 z k , j cos mn cos nny. 

2 — H 



(B) 



ljO_ 



i=k-l 

Z* 

i=l 



j=i-l i-k-1 
.Z Z z. 



iJ. 



mTT x . 
COS 1 



cos nil 



j ft l i =1 



i » J 



mn x . nn y . 

cos i cos J ,i 



(C) 



(D) 



Similarly, Equations (4-6) become: 



^m , n lk 



No (A) term in that on the corner 

-» o 



points, gin mnx. 



(A 



( 10 ) 



No B term due to the fact that 
x Q and x^ are equal to L and sine 

term goes to zero 



+ *7 



i~k-l 

^ AjSlIiA sin mnx 

i=l 2 sin -r- 1 



COS 



nTT 



j-Jl-1 i=k-l 

z 

ri i"i 



mn x . „ nn y . 

v z. . sm l cos J .1 

^ 1 1 j T 



(C) + 
(D)* 



H 



( 11 ) 
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No (A) term because sine terms go to 
zero. 



E ' Z o.i 4z k..1 cos mTT sir) ^ 
J-l 2 H 



(B) 

4 



No (C) term because sine terms go to 
zero . 



j=Jt-l 

£ 

j-l* 



i-k-1 

L 

i=l • 




cos 



mux . 
i 

L 



sin 



nriy 
■ tJ, 

H 



(D) 



No (A), (B) , or (C) terms because sine 
terms to to zero. 



j-J?-l 

£ 

n 



i=k-l 

E 

i=l 



(D) 



z. . sin mn *i 
i » J t. 



sin 



nny 

- j.L 

H 



\ 



( 12 ) 



(. 13 ) 



-s. 



APPENDIX C 

COMPUTER PROGRAM TO DETERMINE COEFFICIENTS FOR SINGLE 
FOURIER SERIES (ONE DIMENSION) FIT OF ORTHOGONAL FUNCTIONS 

Purpose 

The purpose of this program is to determine the Fourier 
coefficients a n and b n , when a function, f(x), is fitted by 
a single Or one-dimensional Fourier series. The secondary 
purpose is to calculate an estimate of the power spectrum 
of the function, which in turn is a characterization of the 
function. 

Language 



Fortran 


IV (IBM 


7040 Computer) 




Symbolic Dictionary 


Variable 


S/A* 


I/O** 


Description 


X 


A 


I/O 


Independent variable data 
point . 


Z 


A 


I/O 


Dependent variable data 
point. 


MZ 


A 


I 


Fixed point value for Z. 


A 


A 


0 


Fourier coefficient. 


■ *• 
B 


A 


0 


Fourier coefficient. 


* S--Single 


Variable 


j A--Array of Variables 


** I--Input 


Variable,; 


0--0utput 


Variable 
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Variable S/A I/O 

C A' 0 

P A 0 

\ 

RAO 



CC S 



Description 



Estimate of the function's 
power spectrum based upon the 
equation : 



C 



2 



2 2 
A ^ + Es 



n n n • 

Estimate of input function 
Z at each data point, created 
by using calculated 
coefficients. • ' 



Residual or difference between 
input function Z and estimated 
function P, at each data point. 

A(N) 2 + B(N) 2 



K S 

EL S 

KOW S 

K00 S 

PIE S 

KKK S 

ZSUM S 

SUMCS S 

SUMSIN S 



I Number of data points. 

I Half-length of sample. 

I Number of samples. 

I Sample index. 

I Mathematical expression input 

as being equal to 3*i4l592o5. 

I Number of harmonics to which 

the Fourier series is to be 
expanded . 

Sum of Z data point values. 

Sum of Sine and Cosine terms.' 

Sum of Sine terms in equation 
to determine coefficients. 



SUMCOS S 

AZERO S 

I S 

PM AX S 



Sum of Cosine terms in equation 
to determine coefficients. 

0 First term of Fourier expansion. 

0 Data point index. 

0 Maximum value of variable P. 



S 



0 Minimum value of variable P. 



PMIN 



Variable S/A I/O Description 



ZMAX 


S 


0 


Maximum value 


of 


variable 


Z. 


ZMIN 


S 


0 


Minimum value 


of 


variable 


Z. 


RMAX 


S 


0 


Maximum value 


of 


variable 


R. 


RMIN 


‘S 


0 


Minimum value 


of 


variable 


R. 


N 


S 


, , 


Harmonic index 


; . 


• 





Program Routine 

The program utilizes data points equally spaced along 
the function from a set origin, and by using Fortran equi- 
valents of Equations- (4) and (6) from Appendix A, generates 
a Fourier series estimate of the function. The Fourier 
coefficients generated are presented as characterizing the 

function, and by combining them according to the "Lee" 

» 

estimate, C(-N)^ ■ (A(N) 2 + B(Nb^)/4, a characterizing power 
spectrum which is independent of phase angle, is' determined. 
When more than 'one sample is run, the C(N) values are grouped 
statistically for analysis. The program checks the values 
of the coefficients by creating an estimate of the original 
function, using Equation (1) from Appendix A. A plot of the 
original function, t^he created estimate, and the residuals, 
or differences between the two is presented. 



n n n n n n r> 
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PROGRAM -FOURIER- AUGUST 1965 

SINGLE FOURIER SERIES PROGRAM (ONE DIMENSION) __ 

PROGRAMMED BY WD ALDENDERFFR 

ONE DIMENSION FIT OF 3MM BEADS DATA. 

400 POINTS 50 HARMONICS 

DIMENSION X (80iT»Z(801 ) * MzTs oTTVa ( 3 0 ) * B ( 30~ ) * c (To - * ' 3 "o 77 pT 8 oT ) » 

1 1 DSH ( 32 ) * TITLE(13)*LO(30)*LC(30) *LR ( 30 ) ♦ R ( 80 1 ) 

CALL CLOCK (0) 

C 

READ (5*102) (TITLE! I )* 1=1*13) 

102 FORMAT ( 13A6) . 

READ (5*1) K* EL 
1 FORMAT ( I10*10X*F10.0) - 

WR I TE ( 6 » 6000 ) 

6000 FORMAT ( l'Hl *20X *2 1HPOWER SPECTRUM VALUES ) 

KOW = 16 

C KOW I S THE NUMBER OF SAMPLES BEING TE STED* 

KOO = 1 ' 

PIE = 3.14159265 
AK = K 

MK = K-l __ • 

AMK = MK 

KKK =50 _ 

C KKK IS THE NUMBER OF HARMONICS* 

DELX = P I E / ( ( AK-1. )/2. ) 

555 X ( 1 ) = -(PIE) 

ZSUM = 0.0 
SUMCS = 0.0 

Z ( 1 ) = 0.0 

DO 203 III = 2 ♦ K 
Z ( 1 1 I ) = 0.0 

X(III) = X(III-f) + DEL X 

203 CONTINUE 

READ( 5*204) ( MZ ( I ) ♦ I =1 ♦ K ) 

204 FORMAT (80 1 1/80 1 1/80 1 1/801 1/8011 /1 1 ) 

DO 205 JJ = 1 » K 

IF (MZ(JJ) - 0 ) 800*800 *801 

800 Z(JJ) = -1.0 
GO TO 205 

801 Z(JJ) = MZ(JJ) 

205 CONTINUE 

DO 12 II = 2 *MK 
ZSUM = ZSUM + Z( II ) 

12 CONTINUE 

SUME =(Z(1) + Z(K))/2. 

AZERO = (2. /AMK) *(ZSUM + SUME) 

DO 4 N = 1 *KKK 
SUMS IN = 0.0 
SOMF = 0.0 
SUMCOS =0.0 
FN = N 
SUM= 0.0 

C ALGEBRAIC INDEX EQUALS FORTRAN INDEX LESS ONE. 

DO 5 I = 2*MK 

SUMCOS = SUMCOS + Z ( I ) *COS ( ( FN* PIE * X ( I ) ) /EL ) 



nnnnn nnnnnn 
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SUMSIN = SUMS I N + Z ( I ) *S I N ( I FN* PIE * XIIJJ/EL) 
5 CONTINUE 

SOME = ( ( Z ( 1 ) + Z ( < ) ) /2 • 0 ) * COS ( FN*P I E ) 



COMPUTING THE VALUE OF THE POWER SPECTRUM ESTIMATE • 
ESTIMATE = A ( N ) SQUARED + R ( N ) SQ UARED DIV IDED BY F OUR 



A ( N ) = 2./AMK * (SUMCOS + SOME) 

BIN) = 2.0/AMK * SUMSIN 

CC = AIN)* AIN) + BIN) * BIN) 

C ( KOO ♦ N ) = CC/4.0 
4 CONTINUE 

WRITE I 6, 687) (C(KOO*N) *N=1»50) 

687 FORMAT! I 1X,10F10.4) ) 

W9 I TE ( 6 * 6005 ) 

6005 FORMAT 1 1H0»11HNEXT SAMPLE) 

GO TO 686 
2 NMAX = (K-l)/2 
DO 8 I = 1*K 
SUMCS = AZERO / 2 , 

DO 9 N = 1 *KKK 
FN = N 

SUMCS = SUMCS + AIN) * COSKFN* PIE * X{I))/ED + 
1BIN) * SINKFN * PIE * XI I )')/EL > 

9 CONTINUE 
PII) = 'SUMCS 
RII) = PII) - Z I I ) 



8 CONTINUE 
GO TO 213 




PLOTTING ROUTINE 



214 ISTAR = -13702990896 
I EYE = 27661634608 

I BLANK = -17997958192 
I DASH = -818089008 
DO 10 I = 1»32 
10 IDSHI I ) = IDASH 

WRITE (6*103 ) (TITLE! I ) ♦ 1 = 1 ♦IS)* 

103 FORMAT I 1H1 , 1 3A6// ♦ 17X * 8HOBSERVED »32X ♦ 10HCALCULATED, 
132X,8HRESIDUAL ) 

WRITE I 6 ♦ 1 04 ) I IDSHI I ) ♦ 1 = 1 .32 ) » ( iDSHTl ) * I «r*32) * ' 

II IDSHI I ) *1=1*32) 

104 FORMAT I 5X*32A1*9X*32A1 *9X*32A1) 

ZM AX = Z (1 ) 

ZM I N = Z I 1 ) 

PM AX = PII) 

PM IN = PII) 

RM AX = R (1 ) 

RM IN = RID 

DO 20 I = 2* K \ 

IF (Z(I).GT. ZMAX) ZMAX = 2(1) 

IF 12(1) .LT. ZMIN) ZM I N = 2(1) 



r>nnnnnnr* nor* 



- 7 - 8 . 



IF (P(I).GT. PMAX ) PMA X = P{I) _ 

IF (P(I).LT. PM I N ) PM I N = P{I) 

IF (R(I).GT. RMAX ) RMAX = R(I)_ 

20 IF (R(I).LT. RMIN) RMIN = R(I) 

WRITE (6,902 ) ' ZMAX »ZM I N »PMAX »PM I N »RMAX »RM I N 

902 FORMAT (IX, AHCHE2 ♦ 6F1 0 • 4 ) 

DO 2 50 I = 1 ,30 

LO ( I ) = IBLANK 
LC ( I ) =i IBLANK 
250 LR( I ) = IBLANK 
KHAF = (K+l ) /2 
DO 30 1= 1 , K 

JZ =( ( ( Z ( I J-ZMIN) / (ZMAX-ZMIN) ) * 29.) + 1. 

JO =( ( ( P ( I ) -PMI N ) / ( PMAX -PM IN))* 29. ) + 1. 

JR =( ( (R( I )-RMIN)/(RMaX-RMIN) ) * 29.) ♦ 1. 

LO (JZ) = I S T A R 
LC (JO) = I ST AP 
LR (JR) = I ST AR 
KK = I - KHAF 

WRITE (6, 105) KK, (LO( JJ) ,JJ = 1 ,30) , KK,(LC( JJ) 73j=l»30) , 
IKK , ( LR ( J J ) » JJ = 1 ,30 ) 

105 FORMAT (1X,IA, 1HI* 30A1V1HI, 5X , I A ♦ 1H I", 30 AT ,TH T, 5X , 

1 1 A , 1 H I »30A1 ,1HI ) 

LO ( JZ ) = IBLANK _ r 

LC ( JO ) = IBLANK 
LR ( JR ) = IBLANK 
30 CONTINUE 



END OF PLOTTING ROUTINE 



213 DO 699 I = 1 , K 

WRITE (6, 698) I , Z ( I ) * P ( I > *R f I ) 

698 FORMAT ( 10X, I5,3F20.8 )_ 

699 CONTINUE 

686 KOO = KOO + 1 

IF (KOO.GT.KOW ) GO v TO A00 
GO TO 555 
A00 CALL CLOCK ( 1 ) 



GROUPING OF C ( N ) FACTORS (POWER SPECTRUM) FOR FREQUENCY 
STUDY. POWER SPECTRUM VALUE = A*A + B*B DIVIDED BY 
A • 0 AS DEFINED IN THE TEXT BY LEE. THIS IS ONLY AN 
ESTIMATE OF THE TRUE SPECTRUM. 



WR I T E ( 6 , 577 ) 

577 FORMAT ( 1H1 , 37HFREQUENCY DISTRl’BUT ION OF“C(N) VALUES) 

WR ITE ( 6 , 575 ) 

575 FORMAT(1HO*3X*1HN»6X»AHINT1»6X, AH I NT 2 , 6X » AHI NT3 » 6X ♦ 
1AHINTA,6X,AHINT5,6X,AHINT6,6X,AHINT7,6X,AHINT8,6X,AHINT9, 
26X ,5HINT10,5X,5HINT11 ,5X,5HINT12) 

WR ITE ( 6 , 576 ) 

576 FORMAT(1HO,9X,AOHLT .001 .001-.01 '"".01-.02 .'02-V03 " » 

150H.03-.0A .0A-.05 .05-. 06 .06-. 07 .07-. 08 ., 

225H08-.09 .09-. 1 GT.1)“ 
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DO 500 N = 1 »KK K 

TOT = 0.0 
TOTS = 0.0 
INTI = 0 
INT? = 0 
INT3 = 0 

INT4 = 0 

INT5 = 0 
INT6 = 0 
INT7 = 0 
INT8 = 0 
I N T9 = 0 

INT10 = 0 ' 

INTI 1 = 0 
INT12 =0 

DO 501 KOO = 1 ♦ KOW 
COEFF = C { KOO ♦ N ) 

TOT = TOT + COEFF 

TOTS = TOTS + <COEFF*COEFF) 

IF (COEFF.GT. . 00099995) GO TO 505 
INTI = INTI + 1 
GO TO 501 

505 IF (COEFF.GT.. 00995) GO TO 506 
INT2 = I NT2 + 1 

GO TO 501 

506 IF (COEFF. GTV.01999'5.) GO TO 507 
I N T 3 = I NT3 + 1 

GO TO 501 

507 I F ( C OE F F . GT . . 02 99 95 ) 00__ JO 5 08 _ 

I NT4 = INT4 + 1 " 

GO TO 501 

508 IF (COEFF.GT. . 039995) GO TO 509 

___INT5 = _I NT 5 +__1 

GO TO 501 

509 IF (COEFF.GT. .049995J__ GO TO__51_0 
INT6 = I NT6 + 1 

GO TO 501 

510 IF(COEFF.GT.. 059995) GO TO 511 
INT7 = I NT 7 +1 

GO TO 501 

511 IF(COEFF.GT. .069995 )_ GO TO 512 
INT 8 = INT8 + 1 

GO TO 501 

512 IF(COEFF.GT. .079995) GO TO 513 
INT9 = I NT9 +1 

GO TO 501 

513 IF(COEFF.GT.. 089995) GO TO 514 
INT10 = INT10 + 1 

GO TO 501 _ _ _ 

514 I F ( COEFF .GT • • 1 ) GO TO 515 
INT11 = INTI 1 + 1 

GO TO 501 " * 

515 INT12 = INTI 2 + 1 
501 CONTINUE 

AKOW = KOW 
AVE = TOT/AKOW 



r\ n r> 
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AVF2 = TOTS/AXOW , ; 

AVE3 = AVE*AVE 

SDEV = SORT ( ( AVE2-AVE3 ) ) 

WRITE ( 6,502 ) N ♦ I NT 1 ♦ I NT 2 * I NT 3 * I NT4 ♦ INT5 ♦ INT6* INT7* I NTS* 

1INT9*INT10*INT11*INT12 

502 FORMAT ( 1H0* 14*12 110) 

WR ITE I 6 *602 ) AVE* SDEV 

602 FORMAT ( 1H0*7HMEAN = ♦ F2 0 . 8 ♦ 1 OX ♦ 1 6HSTANDARD DEV. = *F20.8) 



END OF GROUPING ROUTINE 



500 CONTINUE 

CALL_CLOCK(2) 
CALL EXIT 
END 

SENTRY 
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(*) 
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Y_. 

| r:aoi mz(i) 

/T~^\ 

^ i ,y 








83 



0 



90MK - ((2(1) * 2(K))/2 



(JOS(FN ♦ 



PIS) 



2 




0 
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SIN ILB roURISR 3SRIS3 » L0rTIH3 ROUTINB 




-pj ZHAX * 2(1) 



-> ZMIN *• Z(I) 



> 






*lprfIN 








-J RMAX 




] 


1 


^RMIN 


- R(I) 



L0( I ) - I BLANK* 
LO(I) ' I BLANK 
LR( I ) « IBLANK_ 




250 







KHAP ■=. K f 1/2 







>0 \ 

I * 1 »K , 



JZ - (((Z(I)-ZMIN)/ 

( zmax-zmin) ) • 29.) +• 1. 

JO * (((P(I)-PMIN)/ 

(PMAX-PMIN) ) * 29.)+ 1. 

JR - (((R(I)-RMII®)/ 

(RMAX-RMIN) ) * 29.) + 1. 



LO(JZ) - I3TAR 
LO(JO) - I5TAR 
LR(JR) i I3TAR 



© 




ft 
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COMPUTER PROGRAM TO DETERMINE COEFFICIENTS 
FOR DOUBLE FOURIER SERIES (TWO DIMENSION) 
APPROXIMATION OF A SURFACE AREA 



Purpose 

The purpose of this program is to approximate a surface 
area, flx,y), using a double Fourier series, and to present 
the Fourier coefficients which are a characterization of the 
surface . 



Language 

Fortran IV (IBM yoUO Computer) 



Variable 



X 

Y 

Z 

EL 

H 

KK 



Symbolic Dictionary 



S/A* 

A 



A 

S 

S 

S 



I/O** 

I 



I 

I 



Description 

Independent variable data 
point along x-axis. 

Independent variable data 
point along y-axis. 

Dependent variable data point 

One-half of fundamental per- 
iod along the x-axis. 

One-half of fundamental per- 
iod along the y-axis. 

Number of data points in 
x-direction. 



* S--Single Variable; A--Array of Variables 

** I--Input Variable; 0--0utput Variable 
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Variable ' S/A I/O Description 



LL 


S 


\ 


I 


Number of data points in 
y-direction . 


MMAX 

and 

NMAX 


S 




I 


Maximum number of harmonics 
to which Fourier series is to 
be expanded. Should equal 
one another. 


M 


S 




0 


Harmonic index. 


N 


S 




0 


Harmonic index. 


PIE 


s ' 




I 


Mathematical expression input 
as being equal to 3*1^159265. 


A 


A 




0 


Fourier coefficient. 


B 


A 




0 


Fourier coefficient. 


C 


A 




0 


Fourier coefficient. 


D 


A 




0 


Fourier coefficient. 



Program Routine 

This program utilizes data points equally spaced upon 
a grid, placed upon the surface to be analyzed. An-approxi- . 
mation of the surface is created by expanding a double 
Fourier series of the type discussed in Appendix B. The 
Fourier coefficients are presented as characterizations of 
the surface. 

Input Generator 

Immediately following the listing of the double Fourier 
series program, is a listing of a simple data generator for 
the program. This program uses a random set of coefficients, 
and Equation (2) of Appendix B, to create, in the form of 



90 



punched output, data points simulating a surface area, which 
are the appropriate input if or the double Fourier program. 



no nnn n n on nnnnnnn 



91 



PROGRAM 2D FOUR I AUGUST 1965 

DOUBLE FOURIER SERIES (TWO DIMENSION) PROGRAM 
PROGRAMMED BY WD ALDENDERFER 



DIMENSION X ( 2 5 ) * Y ( 3 1 ) » Z ( 2 5 * 3 1 > * A ( 5 *5 ) 

DIMENSION B (' 5 » 5 ) ♦ C ( 5 ♦ 5 ) » D ( 5 » 5 ) 

CALL CLOCK! 0 ) 

WR I TE ( 6 » 600 ) 

600 FORMAT (lHlVl9H2D FOURIER ANALYSIS ) 

WRITE! 6,623 ) 

623 FORMAT ( 1H0»20X,26HDISPLAY OF CALCULATED 2-D » 
140HC0EFFICIENTS WHICH CHARACTERIZE FUNCTION) 
WRITE(6»654) 

654 FORMAT ( 1 HO »4X ♦ 1HM » 4X ♦ 1 HN » 7X , 6HA ( M »N ) * 14X » 6HB ( M »N ) ♦ 

1 14X » 6HC ( M ♦ N ) * 1 4X » 6HD ( M ♦ N ) ) ~ 

READ ( 5 » 2 ) EL»H»KK»LL 
2 FORMAT (2F20. 8 ♦ 215) 

EL = FUNDAMENTAL LENGTH IN THE X DIRECTION 

H = FUNDAMENTAL 'LENGTH IN THE Y DIRECTION 

KK = NO. OF DATA ^OINTS IN X DIRECTION. 

LL = NO. OF DATA POINTS I N Y D I ft EC TTON . 

READ! 5 » 1 ) ( (X! I )*Y(J) *Z( I*J) * J=1 »LL ) ♦ 1 = 1 »KK ) 

1 FORMAT ( 3F20. 8 ) 

X = INDEPENDENT VARIABLE. 

Y =' INDEPENDENT VAR I ABLE. 

Z = DEPENDENT VARIABLE. 

MMAX = 5 
NMAX = 5 

NMAX AND MMAX ARE THE' NUMBER OF HARMONICS TO WH'ICH“THE 
SERIES IS TO BE EXPANDED. 

PIE = 3.14159265 

AKK = KK-1 

ALL = LL-1 
M<K = AKK 
mll = ALL 
AMK = AKK*ALL 
DO 22 M = l.MMAX 
FM = M— 1 

DO 2 2 N = l" V NMAX 

FN = N-l __ • 

ABLE = 0. 

ABAKE = 0. 

ACHAR = 0.0 
ADELT = 0.0 
BCHAR = 0.0 
BDFLT = 0.0 
CBAKE = 0.0 
CDELT = 0.0 
DDELT = 0.0 
PIM = PIE * FM 
PIN = PIE*FN 



non r> n n nn 



CALCULATE A(0»0> 

ABLE = ( ( ( ( Z ( 1 * 1 ) + Z ( K K * 1 ) ) /2 » 0 ) + ((Z(l.LL) + Z ( KK ♦ LL ) ) /2 .0 ) ) 

1/2.0) * cos(piml«lcos(pinj 

DO 20 I = 2 »MKK 

DO 20 J = 2 »MLL _ 

ARGX = PIM * ( X ( I) /EL) 

ARGY = PIN * ( Y ( J ) /H) 

CX = COS (ARGX) 

CY = COS ( ARGY ) : 

SX = SIN(ARGX) 

SY = SIN (ARGY) 

I F ( I.GT.2) GO TO 400 

CALCULATE COEFFICIENTS BY PARTS. 

ABAKE = A B A K E + UT iff* J f T Z ( KK . JTT72T6 ) ~ *CX *CY~ 

CBAKE = CBAKE + ((Z(1*J) + Z ( KK ♦ J ) ) /2 • 0 ) * CX * SY 

400 IFU.GT.2) GO TO 401 

ACHAR = • ACHAR+ ((Z(I»1) + Z ( I *LL ) ) /2.0) *CX* CY 
BCHAR = BCHAR + ( ( ZVl Vi ) + Z ( I * LL ) ) /2.0 )*SX*CY~ 

401 ADELT = ADELT + Z(I*J) *CX*CY 

BDELT - BDELT + Z ( I * j ) *SX*CY 

CDELT = CDELT + Z(I*J) * CX * SY 
DDELT = DDELT + Z(I.J) * SX * SY 
20 CONTINUE 

A(M*N) = ( ABLE+ ABAKE+ ACHAR + ADELT) * ( 4 .0/ ( AKK*ALL ) ) 

B(M*N) = (BCHAR + BDELT) M4.0/AMK) 

C ( M * N ) = (CBAKE + CDELT) * (4.0/AMK) 

D(M.N) = (DDELT) * (4.0/AMK) ; 

22 CONTINUE 



WRITE VALUES OF COEFFICIENTS. 



DO 25 M = ltMMAX 
MM = M-l 

DO 25 N = 1 VNMAX 
NN = N-l 

WRITE (6 .27") MM»NN»A(M»N) »B(MVN) »C(MtN) »D(M*N) 

27 FORMAT ( 1 HO »2 I 5 »4F20.4 ) 

25 CONTINUE 

_ CALL CLOCK ( 1 ) _ _ _ _ 

CALL EXIT 
END _ 

SENTRY 



nnnnnnnnn 
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PROGRAM 2D GENERATOR 

program to generate dependent variable z as a function 

OF TWO INDEPENDENT VARIABLES X AND Y. PROGRAM 
USES DOUBLE FOURIER EQUATION AND GIVEN COEFFICIENTS*. 
COEFFICIENTS USED IN THIS PARTICULAR PROGRAM ARE 
TAKEN FROM PRESTON-HARBAUGH PAPER. 

PROGRAMMED SEPTEMBER 1965. 

DIMENSION A(5*5)*B(5*5) ♦ C ( 5 » 5 ) * D ( 5 ♦ 5 ) ♦ T (25.31 ) 

CALL CLOCK (0) 

DO 100 N = 1 * 5 
DO 100 M e l ♦ 5 
FN=N 

fm=m 

A(M.N) = (.3 + • 02*FN + . 1 1*FM ) *1 00 • 

B(M.N) = (.2 + .04*FN + .08*FM)*100. 

C ( M « N) = (.1 + ,06*FN + .05*FM ) *1 00 • 

D(M.N) = (.08 + .05*FN + .02*FM) *100. 

100 CONTINUE 



DO 


101 


M 


*i 


♦ 5 


C ( M 


»1 ) 


3 


0. 


0 


B ( 1 


*M) 


s 


0. 


0 


D ( 1 


♦ M ) 


r 


0. 


0 


101 D ( M 


♦ 1 ) 


= 


0. 


0 



PI = 3 . 14159265 
SUMZ = 0.0 
AH = 15. 

AL = 12. 

WR I TE ( 6 ♦ 984 ) 

984 FORMAT ( 1H1 »29H INPUT VALUES OF COEFFICIENTS.) 

WRITE( 6.982 ) 

982 FORMAT ( 1H0.5X*81HM N A(M.N) B ( M » N) 

1 C(M.N) D ( M »N ) ) 

DO 300 N = 1*5 

DO 300 M = 1*5 

WRITE(6*981)M*N*A(M.N).B(M*N)*C(M*NJ*D(M*N) 

981 FORMAT (2X.2I5.4F20.8) 

300 CONTINUE 

PIL = PI /AL ' ' 

P I H=P I /AH 
XC = -12. 

YC = -15. 

XI NC =, 1.0 

YINC = 1 .0 ______ __ 

WR ITE ( 6 ♦ 983 ) 

983 FORMAT ( 1H1 .11X»24HX( I ) Y£J) * 

122H Z ( I ♦ J ) * 1 3H I J) 

DO 410 I = 1*25 
DO 405 J * 1*31 
DO 200 N=1 ♦ 5 
X = XC 
Y = YC 
FNT = N-1 

PINH * P IH*FNT 



DO 200' M =1 ♦ 5 
fmt=m-i 

P I ML = PIL*FMT 
IF (N-l ) 1 40 ♦ 1 20 * 1 40 

120 IF (M-l ) 15.0*1.30*150^ 

130 AMDA=.25 
GO TO 170 

140 IF (M-l ) 160*150*160' 

150 AM DA = • 50 
GO TO 170 

160 AMDA=1.00 : 

170 ARGX = P IML*X 

ARGY = P I NH*Y 

CX = COG (ARGX) 

CY = COS (ARGY) .. ___ 1 

SX = SIN(ARGX) 

SY = SIN (ARGY) _ . 

ACC = A (M»N)*CX*CY 

DSC = B ( M * N ) *SX*CY _ _ _ ; „ 

CCS = C ( M ♦ N ) *CX*SY 
DSS = D ( M ♦ N ) *SX*SY 

TS = AMDA* ( ACC + BSC + CCS + DSST " 

_ZC = ZC + TS 

200 CONTINUE 

WRITE (6, 980) X * Y * Z C » I ♦ J 
980 ' FORMAT ( 1 X * 3F20 . 8 *215)” 

WR ITE ( 7 * 950 ) X*Y*ZC 

950 FORMAT (3F20. 8 ) \ 

L T( I , J) =_ZC 

ZC = 0.0 

405 YC = YC + YINC _ . 

YC = -15. 

410 XC = XC + XINC 

t ( 1 * 1 ) = ( T ( 1 ♦ 1 ) + T( 1 *31 ) + T ( 25 * 1 ) +"T ('25*31 ) )/4.0 

T( 1 *31 ) = 0 .0 ^ ■_ 

" T(25*l ) = 0.0 
T ( 25 *31 ) = 0.0 
DO 700 J = 2*30 

T ( 1 ♦ J ) = ( T ( 1 ♦ J ) _ +_ T_( 2 5 * J ) ) / 2 * 0 

T ( 25 ♦ J ) = 0.0 

_700_ CONTINUE 

DO 701 I = 2*24 

T ( I ♦ 1 ) = (T( I *1 ) +_ J_( 1*31 ) )/2.0 _ 

T( I *31 ) = 0.0 
... 701 CONTINUE 

DO 702 I = 1*25 

DO 702 J = 1*31 

SUMZ = SUMZ + T( I * J) 

702 CONTINUE 

A00 = (4.0*SUMZ )/720. 

WR I T E ( 6 » 9999 ) SUMZ* AOO 

9999 FORMAT (1H1*18HSUM OF Z VALUES '♦ F20V8V5 X * 9HA ( 0 * 0 ) V VF20.8 ) 

CALL CLOCK (1) 

490 CALL EXIT 
END 

.. TENTRY 



DOUBLE FO RISp ‘JKRIK3 PROGRAM 
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j 3T A RT j 



T VRITfil nrLE 

^ and HEADINGS 



^READi EL, H, KK, uT| 

j 

^RtAD « X(I), Y(J), Z(I,J)^ 



NMAX = NO. 


HARHOMoi 


MMAX ' NO. 


HARMON I 03 


"1 


r ■■ 


^PIB = 5.14I59265J 




BOHAR * 0. 
| BDELT s. 0. 
* OBAKE x 0. 

I 

ODEIT - 0. 






(& 



% 
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l!) 



aoslt = adelt ♦z(i.j) * ♦ oY 

-OELT = bOSLT ♦- Z( 1 , J ) * 3X • OY 

ODELT - OOELT + 7,{ I , J ) * OX * 3Y 

D05LT * DOSLT^-Z(I.J) * 3X * 3Y 



Q 
1 ® 



— J I 
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COMPUTER PROGRAM TO DETERMINE COVARIANCE 
FUNCTION AND POWER SPECTRUM 



Purpose 

The purpose of this program is to determine the 

/ * 

covariance function ^nd the power spectrum of a one- 
dimensional function, f(x). The secondary purpose is to 
plot the power spectrum and the confidence bands of the 
spectrum, to facilitate analysis of the function. 



• » 

Language 

Fortran IV (IBM 70^0 Computer) 



Variable 



Symbolic Dictionary 



S/A* I/O** 



Description 



X 



I 



A' I 



S 



\ 



Value of the function. at each 
data point. Read in as fixed 
point, mx. 

Series index. 



K 

C 

L 

U 



S 

A 

A 



0 

0 



Lag Index. 

Covariance estimate. 

Raw estimate of the power 
spectrum. * 

Smoothed estimate of the 
power spectrum. 






* S--Single Variable; A--Array of Variables 

** I--Input Variable: 0--0utput Variables 



Variable 


S/A 


I/O 


Description 


N 


S' 


I 


Number of data points. 


M 


s 


I 


Maximum lag. (Less than N/3) 


EDF 


A 


I 


Equivalent degrees of freedom. 


TL95 


A 


I 


Table of confidence limits, 
lower 95 per cent.. * 


TL90 


A 


I 


Table of confidence limits, 
lower 90 per cent. 


TU95 


A 


I 


Table of confidence limits, 
upper 95 per cent. 


TU90 


A 


I 


Table of confidence limits, 
upper 90 per cent. 


KOW 


S 


I 


Number of samples. 


KOUNT 


S 


-- 


Sample index. 


LA 


A 


I 


Storage for data identifi- 
cation . 


LAG 


S 


0 


Degree of offset by which 
the function is correlated 
with itself, to determine 
covariance estimate. 


UMAX 


S 


— 


Maximum value of smoothed 
power spectrum in each fre- 
quency band. 


UMIN 


S 


-- 


Minimum value of power 
spectrum . 


SPANA 


S 


-- 


90 per cent confidence band 
of smoothed power spectrum. 


SPANB 


S 


— 


95 per cent confidence band. 



Program Routine 

Thi's program utilizes data points which represent the 

value of the function to be analyzed at equally spaced 

« 

intervals along the function. Using a numerical integration 
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techniaue, the trapezoidal rule, the values of the 

covariance estimate are determined according to Equation (3*25) 

of this paper. A rough estimate of the power spectrum is 

computed using Equation ( 3 . 26 ) ard smoothed according to 

Equation (3*2?). The table of confidence limits presented by 
77 

Granger is read into the' program and the ninety and 
ninety-five per cent confidence bands are determined. 

The smoothed power spectrum with these confidence bands is 
displayed graphically. 



UUUUULMJUUUUUUUUUUUUUUUUUU , i , uuu 
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PROGRAM COVAR FOR COMPUTING COVARIANCE FUNCTION 
AND POWER SPECTRUM. PROGRAMMED SEPT. 1965 BY 
FW PRESTON AND WD ALDENDERFER. 



NOMENCLATURE 

< = lag index 

I = SERIES INDEX 

J = WEIGHTING FACTOR INDEX 

C ( K ) = COVARIANCE ESTIMATE 

L ( J ) = RAW ESTIMATE OF POWER SPECTRUM 

U(J) = SMOOTHED ESTIMATE OF POWER SPECTRUM 
M = MAXIMUM LAG (LESS THAN N/3) 

N = NUMBER OF DATA POINTS 

LA ( I ) = STORAGE FOR DATA IDENTIFICATION 

EDF(I) = EQUIVALENT DEGREES OF FREEDOM 

TL95 ( I ) = TABLE OF C ONFIDENC LIMI TS ~ LO WER- 95PERCENT_ 

TL90(I)= - LOWER- 90PERCENT 

TU90 ( I ) = - UPPER 90 PERCENT 

TU95 ( I ) = - UPPER 95 PERCENT 

KOW IS THE NUMBER OF SAMPLES 



D I MENS ION X ( 3200 ) » C ( 30 1 ) » L ( 302 ) » U( 301 ) » EDF (13) * TL95 ( 13 ) , 
1TU90( 13 ) »TU95 ( 13 ) , LA ( 13 ) ,J0W1 (53 ) ,JOW2 (53) »TL90( 13 ) .MX (3200 ) 
REAL L 

CALL CLOCK (0) 

KOW = 2 
KOUNT = 1 

READ (5*110) (EDF ( I ) ,TL95( I' ) » TL90 ( I ) *TU90< I ) *TU95H ) *1 = 1*13 )‘ 
110 FORMAT ( F5. 1 »4F 5.. 2_)_ 

400 READ (5*99) (LA (I) ,1=1*13) 

99 FORMAT (13A6) 

READ (5,100) N ♦ M 

100_FORMAT _( 2 I 5 ) 

READ (5,101) (MX ( I ) ,1 = 1 ♦ N ) 

101 FORMAT (8011 ) ' 

DO 700 11= 1 *N 

X ( 1 1 ) = MX ( 1 1 ) . 

700 CONTINUE 

DO 10 I = 1*N 

IF ( X C I ) - 0. ) 11 *11,10 

11 X ( I ) = -1.0 • * . 

10 CONTINUE 

"COMPUT ING "COVArTaNCE 'ESTIMATE "cTkT" T 



FM = M 
MP 1 = M+l 
FN =N 
J = 0 

AK = 2.0*FN/FM 
DO 910 I = 1,13 
J = J+l 

I F ( AK. LT • EDF ( I ) ) GO TO 320 
910 CONTINUE 



n r, n n no nr»n ann 
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320 J = J-l _ 

CHG = (AK-EDF( J) ) / { EOF ( J + l )-EDF( J) ) 

T90L = T L90 ( J ) + ( TL90 ( J + l ) -TL90 ( J ) ) *CHG 
T95L = TL95 ( J ) + ( TL9 5 ( J+ 1 ) -TL95 ( J ) ) *CHG 
T95U = TU95(J) + ( TU95 { J+l ) -TU95 ( J ) ) *CHG 
T90U = TU90 ( J ) + (TU90( J+1)-TU90( J) )*CHG 

DO 210 KK = 1 ♦MPl 

K = KK- 1 
Sl= 0.0 
NMK = N-K 

FNMK = NMK . 

KP1 = K+l 

DO 205 1=1 ♦ NMK 

IPK = I+K 

205 SI = SI + XII )*X( IPK) 

52 = 0.0 

DO 206 I=KP1.N 

206 S2 = S2 + X( I ) 

53 = 0.0 

DO 207 I=1»NMK 

207 S3 a S3 + Xm 

210 C ( KK ) = (SI - (S2*S3 > /FNMK ) /FNMK 

COMPUTATION OF RAWEST I MATE OF, POWER SPECTRUM L ( J ) 



MN 1 = M-l . • ■ 

DO 230 JJ = 1 *MP 1 
J= JJ-1 
F J= J 
SI =0.0 

DO 220 K = 1 »MN1 

fk=k 

220 SI = SI + C(K+1)*C0S(3.1415926*FK*FJ/FM) 

230 L ( J J+ 1 ) =(2.0*S1 V C(l) +C( M+l ) *C0 S ( 3 . 1 4 159 2 6*FJ) )“/6 . 2831952 
L ( 1.) =L ( 3 ) 

L ( M+3 ) = L(M+1) 

COMPUTATION OF SMOOTHED POWER SPECTRUM ESTIMATE' 



MP2 = M+ 2 
DO 240 J =2»MP2 

240 U ( J ) .= • 25*1 (J-l ) + . 50*L ( J ) WV25*L(J+1) 
“W R I T E~HEAD I NGS 



WR ITE ( 6 » 103 ) 

103 FORMAT ( 1H1 .36X.25HPOWER SPECTRUM ESTIMATION/ 
137X.24HFOR DISCRETE TIME SER I ES/3 IX ♦ 
237HUSING TUKEY HANNING WEIGHTING FACTORS///) 

WRITE{6»104) (LA( I ) ♦1=1.13) 

104 FORMAT ( 1X.13A6/// ) 

WR ITE( 6 ♦ 106 ) 



WRITE POWER SPECTRUM AND TOLERANCE ESTIMATES 



n n n non 



LAG = 1-2 
U90L = T90L*U( I ) 

U95L = T95L*U(I) 

U90U = T 90U*U ( I ) 

U95U = T95U*U( I ) 

250 WRITE(6*105) LAG*C( II ) *L( I) * U90L*U( I ) * U95U »U95 L *U ( I ) *U9QU 
105 FORMAT ( 1 X » I 3 ♦ 8 ( 2 X * F 1 0 • 4 ) ) 

PLOTTING ROUTINE 



764 ISTAR =-13702990896 
IMINUS = 29809118256 
IPLUS = 17997958192 

I BLANK = -17997958192 

II = 27661634608 
WRITE (6*103) 

WRITE (6*1 04) ( L A ( I ) *1=1*13) 

WRITE (6*107) 

107 FORMAT ( 19X*23HSMOOTHED POWER SPECTRUM ♦ 34X ♦ 

123HSMOOTHED POWER SPECTRUM /15X* 

230HWITH 90 PCT. CONFIDENCE LIMITS, 27X* 

330HWITH 95 PCT. CONFIDENCE LIMITS/6X, 

451H0 123 4. 56789 0* 

56X »51H0 123456789 0) 

106 FORMAT ( lH0»39Xr‘ 

150HSMOOTHED SPECTRUM WITH INDICATED CONFIDENCE LIMITS/ 

234X *7H90 PCT . * 17X ♦ 7H90 PCT.*5X,7H95 PCT.V17XV7H95“PCT./8X," 
31 OH CO VAR IANCE*5X*3HRAW*9X* 5H LOWER *19X*5HUPPER»7X* 
45HlOWER*19X*5HUPPER/1X,3HLAG » 5X * 8H5PECTRUM T4X » 8HSPECTRUM ♦” 
56X,5HLIMIT*6X*7HAVERAGE*6X*5HLIMIT* 7X*5HLIMIT *6X ♦ 7HAVERAGE ♦ 
66X , 5HL I M IT ) 

PLOTTING OF SPECTRUM - ' 

DO 280 1 = 1 *53 T“ 

JOW1 ( I ) = IM I NUS 
280 JOW2 ( I ) = I M I NUS 

WRITE (6, 108 ) ( JOW1 ( I ) ,1 = 1*53 ) » ( JOW2 ( I ) ♦ 1 = 1 ,53 ) 

108 FORMAT (5X*53A1 * 5X * 53AT ) 

DO 290 1=2*52 

JOW1 ( I ) = IBLANK 
290 JOW2 ( I )= IBLANK 
JOW1 (1 )=II 
JOW1 ( 53 )=I I 
JOW2 ( 1 ) =11 
JOW2 ( 53 ) =1 1 
’ UMAX = U(2 ) 

UM I N = U ( 2 ) 

DO 300 I = 3 * MP2 
I F ( u ( I ) .GT.UMAX) UMAX=U(I) 

IF(U( I ) .LT.UMIN) UMIN=U ( I ) 

300 CONTINUE 

SPANA = ALOG 1 0 ( UMAX*T 90U ) - ALOGl 0 ( UM I N*T 90L ) 

SPANB = ALOG 10 ( UMAX*T95U ) - ALOGl 0 ( UM I N*T95L ) 
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ZA = ALOG10(UMIN*T90L ) 

ZB = ALOG10(UMIN*T95L) 

DO 310 I = 2 »MP? 

I A 1 = ( ( ALOGlO I U( I ) *T90L ) - ZA ) /SPANA )*50. + 2.0 
I A2= ( ( A LOG 1 0 (U( I ) )-ZA)/SPANA)*50. + 2.0 
I A 3= ( ( ALOGlO (U( I ) *T90U ) - ZA ) /SPANA ) *50 • + 2.0 

JOWl ( I A 1 ) = IMINUS 

J0W1 ( I A2 ) = I STAR 
J0WKIA3) = I PLUS 

I B 1 = ( ( ALOGlO ( U ( I ) *T95L ) - ZB ) /SPANB )*50.+ 2.0 
IB2 = ( (ALOG10(U(I) )-ZR)/SPANA)*50. + 2.0 
I B 3= ( (ALOG10(U( I )*T95U)-ZB ) /SPANB )*50. + 2.0 

J0W2UB1) = IMINUS _ 

J0W2 ( I B2 ) = ISTAR 
J0W2 ( I B3 ) = I PLUS 
LAG =1-2 

WRITE (6. 109) LAG » ( JOWl ( J) , J = 1 * 53 ) ,L AG , ( J0W2 ( J) , J = 1 ,53 ) 
109 FORMAT (IX, I3,1X,53A1,1X, I3»1X,53A1) 

JOWl ( IA1 ) = I BLANK_ 

J0WKIA2) = I PLANK 
JOWl ( I A3 ) = I BLANK 
J0W2 ( IB1 ) = IBLANK 
J0W2 ( I B 2 > = I BLANK 
J0W2 ( I B3 ) ** IBLANK 
310 CONTINUE 
C 

C END OF PLOTTING ROUTINE 



167 IF (KOUNT.EQ.KOW) GO TO 765 



CALL CLOCK (2) 
KOUNT = KOUNT + 1 



GO TO 
765 CALL 
END 

SENTRY 

2. .05 

3. . 12 


400 

EXIT 










.10 2.30 
.20 2.08 










2.99 

2.60 






• 


4. 


.18 


.26 1.94 


2.37 








5. 


.23 


.32 1.85 


2.21 








6. 


.27 


.37 1.77 


2.10 








8. 


.34 


.44 1.68 


1.94 








10. 


.39 ' 


' .49 1.60 


1.83 








12. 


.43 


.53 1.55 


1.75 








15. 


.48 


.57 1.48 


1.66 








20. 


.54 


.62 1.42 


1.51 




4 




30."" 


.62 


.69 1.34 


1 .46 








50. 


.69 


.75 1.26 


1.34 








100.' 


".77 


.82 1.18 


1.22 
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PROOiMM TO D6TBRMIM OOVARIANOE 
PIT30T I OK AND FOvfER SPECTRUM 





START 

X 




KOV * NO. SAMFLE9 






K0UN1 


• - 1 



RBADi Enp(I), TL95(I), TL90(I), 
TU90(I), TU95(I) 




^ *»00 REAOi TITLE 




“1 



_J 



~i 
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\ 



107 







FNKK * KMX 






KPI - Kfl 







205 
I * 1 ,NMK 



/ 



31 » H 1 + X(I) * X(IFK) 




I 



|~ 32 ^ 0 « ] 



A 



206 v 
1 =» kpi, n r 



92 * 32 X( I ) 




207 

i - 1 ,rmk y j 



I 

A- 1 ? 

X I 

( 32 * 3 }) / TOtt / FKMI)| 



o(xx) « (31 - (32*35) / fhmk / ran) 




l 



i 


M -1 







C 2J0 
JJ * 1 ,MP 1 



(2) (°) 
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APPENDIX D 



PICTURES 




Sample, X Three Millimeter Beads 
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Sample Y Two Size (4mm. and 5mm.) Beads 
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lip® 












Photomicrograph of porous medium 
(sandstone) 92X. Courtesy of Humble Oil and Re-, 
fining Company. Porosity 26.3 per cent. The line 
is 0.17 cm in length. 
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POWER SPEC I RUM ESTIMATION 
FOR OISCRETF TIME SERIES 
USING TUKEY-MANNING WEIGHTING FACTORS 



TWO SIZE BE AO$ *0 ROWS MIXEO LAG * 40 

fi) 



swiothep spectrum with inoicateo confidence LIMITS 









90 PCT. 




90 PCT. 


95 PCT. 




95 PCT, 




COVAr I AnC6 


RAW 


LOWER 




UPPER 


LOWER 




UPPER 


LAG 


SPECTRUM 


SPECTRUM 


LIMIT 


AVERAGE 


LIMIT 


LIMIT 


AVERAGE 


LIMIT 


0 


0.9409 


0.7021 


0.6647 


0.7353 


0.7911 


0.6367 


0.7353 


0.7970 


I 


0.6656 


0. 7684 


0.6291 


0.6959 


0.7400 


0.6027 


0.6959 


0.7544 


2 


0.4372 


0.5449 


0.5591 


0.6185 


0.6655 


0.5356 


0.6185 


0.6705 


3 


0.2795 


0.6159 


0.5069 


0.5607 


0.6033 


0.4856 


0.5607 


• 0.6070 


4 


0. 1660 


0.4663 


0.4365 


0.4829 


0.5196 


0.4182 


0.4829 


0.5234 
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0.0853 


0.3831 


0.3751 


0.4149 


0.4464 


0.3593 


0.4149 


0.4497 


6 


0.0401 


0.4270 


0.3815 


0.4220 


0.4541 


0.3654 


0.4220 


0.4574 


7 


0.0053 


0.4507 • 


0.3580 


0.3960 


0.4261 


0.3430 


0.3960 


0.4293 
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-0. 0057 


' 0.2556 


0.2522 


0.2790 


0.3002 


0.2416 
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0.1543 
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0.1707 
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0.1313 


0.1424 
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0.1165 


0.1289 
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0.1116 
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0.1397 
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0.1001 


0.1 195 


0.1286 


0.1035 
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0.1296 
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0.0475 


0.0036 
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0.1052 


0.1132 


0.0911 


0.1052 


0.1140 


14 
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0.1227 
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0. 1105 


0.1109 


0.0957 


0. 1 1C5 


0. 1 190 
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-0.0066 


0. 1130 


0.0919 
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0.1094 


0.0800 


0.1017 


0.1102 
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-0.0170 


0.0500 


0.0652 


0.0721 


0.0776 


0.0624 


0.0721 


0.0782 
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0.0594 


0.0531 


0.0588 


0.0633 


0.0509 


0.0588 


0.0637 


16 
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0.0583 


0.0530 


0.0595 


0.0640 


0.0515 


0.0595 


0.0645 
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0.0515 


0.0569 


0.0613 


0.0493 
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0.0477 


0.0514 
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0.0476 


0.0383 


0.0442 
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0.0436 


0.0351 
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0.0268 
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0.0256 


0.0296 


0.Q321 


20 


0.0032 ' 


0.0235 


0.0250 


0.0276 


0.0297 


0.0239 


0.0276 


0. 0299 


29 


-0.0079 


0.0343 


0.0284 


0.0314 


0.0338 


0.0272 


0.0314 


0. 0340 


30 


-0.0253 


0.0335 


0.0275 


0.0304 


0.0327 


0.0263 


0.0304 


0.0329 


31 


-0.0262 


0.0204 


0.0223 


0.0247 


0.0266 


0.0214 


0.0247 


0.0268 


32 


-0.0133 "" 


0.0246 


0 . 0207 


0.0229 


0.0246 


0.0190 


0.0229 


0.0248 


33 


-0.0067 


0.0221 


0.0201 


0.0222 


0.0239 


0.0192 


0.0222 


0.0241 


34 


-0.0127 


0.0201 


0.0189 


0.020$ 


0.0224 


' 0.0101 


0.0209 


0.0226 


35 


-0.0162 


0.0212 


0.0200 


0.0221 


0.0237 


0.0191 


0.0221 


0.0239 


36 


-0.0273 


0.0258' 


0.0215 


0.0238 


0 .0256 


0.0206 


0.0236 


0.0257 


37 


-0.0333 


0.0222 


0.0221 


0.0244 


0.0263 


0.0211 


0.0244 


0.0265 


36 


-0.0330 


0.0274 


0.022$ 


0.0254 


0.0273 


0.0220 " 


0.0254 


0.0275 


39 


-0.0391 


0.0245 


0.0209 


0.0232 


0*0249 


0.0201 . 


0.0232 


0.0251 


40 


-0.0413 


0 » 0 1 6 3 


0.0184 


0.0204 


0.0219 


‘ 0.0176 


0.0204 


0.0221 



\ 

L 



POWER SPECTRUM ESTIMATION 
FOR 0 1 SC RE T F TIME SERIES 
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POWER SPECTRUM ESTIMATION 
FCR DISCRETF I IME SERIES 
USING TUKEY-HANNING WEIGHTING FACTORS 



TWO SIZE 8EAOS AO ROWS MIXEO LAG * AO 

( 1 ) 



SMOOTHED SPECTRUM WITH INDICATEO CONFIDENCE LIMITS 









90 PCT • 




90 PCI. 


95 PCT. 




95 PCT. 




covariance 


RAW 


LOWER 




UPPER 


LOWER 




UPPER 


LAG 


SPECIRUM 


SPECTRUM 


LIMIT 


AVERAGE 


LIMIT 


LIMIT 


AVERAGE 


LIMIT 


0 


0. 97AR 


0.6A57 


0 ♦ 5 5 A 9 


0.6138 


0.6605 


0.5316 


0.6138 


0.665A 


1 


0*7058 


0.5820 


0.5712 


0.6318 


0.6798 


0 . 5 A 72 


0.6318 


0 • 6 8 A 9 


2 


0 • A 705 


0.7177 


0.5AA6 


0.6025 


0 • 6 A 8 3 


0.5217 


0.6025 


0.6531 


3 


0.2788 


0.3926 


0 • A 995 


0.5526 


0. 59A6 


0.A785 


0.5526 


0.5990 


A 


0.1308 


0.7075 


0.5338 


0.59C5 


0.635A 


0 . 5 1 1 A 


0.5905 


0 • 6 A 0 1 


5 


0.0277 


0.55A5 


0.5319 


0.588A 


0.6331 


0.5095 


0.588A 


0.6378 


6 


-0.0383 


0.5369 


0 • A 693 


0.5191 


0.5585 


0.AA95 


0.5191 


0.5627 


7 


-0.0681 


0 . A A 80 


0.3783 


0 . A 1 8 A 


0. A502 


0 • 3 6 2 A 


0 . A 1 8 A 


0 • A 536 


8 


-0.0766 


0.2A08 


0.2665 


0.29A8 


0.3172 


0.2553 


0.29A8 


0.3195 


9 


-0.0500 


0.2A95 


0. 1885 


0.2085 


0.22A3 


0.1805 


0.2085 


0.2260 


10 


-0.0153 


0.09A1 


0.1307 


0. 1 AA6 


0.1556 


0.1252 


0. 1AA6 


0. 1567 


11 


0.0219 


0. 1 A06 


0.1060 


0.1172 


0.1262 


0.1015 


0.1172 


0.1271 


12 


0.0A35 


0.0937 


0. 1035 


0. 1 1 A 5 


0. 1232 


0.0991 


0.11A5 


0. 12A1 


13 


0.0376 


0. 1299 


0.0986 


0.1090 


0.1173 


0.09AA 


0.1090 


0.1182 


1A 


0.0379 


0.0826 


0.0871 


0.0963 


0.1036 


0 • 08 3 A 


0.0963 


0. 10AA 


15 


0* 02 1 A 


0.0902 


0 • 0 7 A 5 


0.082A 


0.0887 


0 • 0 7 1 A 


0. 082 A 


0. 089 A 


16 


0.01A9 


0.0668 


0.0619 


0.0685 


0.0737 


0.0593 


0.0685 


0.07A2 


17 


0.016A 


0.0501 


0.0501 


0.0555 


0.0597 


0.0A80 


0.0555 


0.0601 


18 


0.005A 


0.05A8 


O.OA83 


0.053A 


0.05 7 A 


0.0A62 


0.053A 


0.0579 


19 


0.0020 


0.0537 


0.0AA7 


0.0A95 


0.0532 


O.OA28 


0.0A95 


0.0536 


20 


-0.0065 


0.0356 


0.0395 


O.OA37 


0.0A70 


0.0378 


O.OA37 


0.0A7A 


21 


-0.01C5 


0.0500 


0.038A 


0. 0A2 A 


O.OA57 


0.0368 


0.0A2A 


0.0A60 


22* 


-0.0158 


0.03A3 


0.0369 


0.0A08 


0.0A39 


0.035A 


0.0A08 


0. OA A 3 


2 3 


-0.0268 


O.0AA8 


0.0355 


0.0393 


O.OA23 


0.03AO 


0.0393 


0.OA26 


2 A 


-0.0227 


0.0333 


0.0322 


0.0356 


0.0383 


0.0308 


0.0356 


0. 0386 


25 


-0.0250 


0.0309 


0.0280 


0.0310 


0.0333 


0.0268 


0.0310 


0.0336 


26 


-0.0252 


0.0289 


0.0281 


0.0311 


0.033A 


0.0269 


0.0311 


0.0337 


27 


-0.0305 


0.0358 


0.0283 


0.0313 


0.0337 


0.0271 


0.0313 


0. 03 AO 


"28 


-0.0A66 


0.0250“ 


0.0275 


0.03CA 


0.0327 


0.0263 


0.030A 


0.0329 


29 


-0.0500 


0.0358 


0.0262 


0.0290 


0.0312 


0.0251 


0.0290 


0.0315 


30 


— 0 • 0 A 60 


0.0196 


0.0219 


0.02A2 


0.0261 


0.0210 


0.02A2 


0.0262 


31 


-0.0388 


0.0219 


0.0180 


0.0199 


0.0215 


0.0173 


0.0199 


0.0216 


32 


-0.0385 


0. 0 1 6 A 


0.020A 


0.0226 


0.02A3 


0.0196 


0.0226 


0.02A5 


33 


— 0.0 1 A 2 


0.0358 


0.0261 


0.0289 


0.0311 


0.0250 


0.0289 


0.0313 


3A 


0.0189 


0.0276 


0 • 026 1 


0.0289 


0.0311 


0.0250 


0.0289 


0.0313 


35 


0.0533 


0.02A6 


0.0223 


0 . 02 A 7 


0.0266 


0 .02 1 A 


0 . 02 A 7 


0.0268 


36 


0.0777 


0.0220 


0.0208 


0.0230 


0.02A7 


0.0199 


0.0230 


0.02A9 


37 


0.0805 


, 0.0232 


0.019A 


0.0215 


0.0231 


0.0186 


0.0215 


0.0233 


"38 


0.0726 


6.0175 


0.0189 


0.0210 


0.0225 


0.0181 ~ 


0.0210 


0.0227 


39 


0.0557 


0.0256 


0.0I9A 


0.0215 


0.0231 


0.0186 


0.0215 


0.0233 


AO 


0.0262 


0.0171 


6.0193 


0.02 1 A 


0.0230 


0.0185 


0.021A 


0.0232 
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TWO SIZE BEADS *0 ROWS MIXED LAG - AO 

M 

SMCOTHEO POWER SPECTRUM SMCOTHEO PCWER SPECTRUM 

WITH 9C PC T • CCNriDENCE LIMlIS WITH 95 PCI. CONFIDENCE LIMITS 
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PCWFR SPECTRUM ESTIMATION 

ftp discrete time series 

USING f l K c Y-MANN I NG WEIGHTING FACT OR s 



THREE Sl/F PE ACS 40 ROWS MIXED l AC * 4 0 

0 ) 



SMOOTHED SPECTRUM WITH 1N01CATE0 CONFIDENCE LIMITS 



_ 






90 PCI. 




90 PCI. 


95 PCI. 




95 PCI. 




COVAR l ANCE 


RAw 


LCwT « 




UPPE R 


LIJwbR 




UPPER 


LAG 


SPEC T RUM 


SPEC TRUM 


LIMIT 


AVERAGE 


l im i r 


LIMIT 


AVERAGE 


LIMIT 


0 


0.9035 


0.4561 


0.4 3 71 


0.4P 35 


0.5203 


0 • 4 1 A 7 
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PCWER SPECTRUM ESTIMATION 
FCR 0 I SC RE T F TIME SERIES 
USING TUK FY-HANN ING WEIGHTING FACTORS 



3MM BEADS 40 RCWS MlXfO LAG = 40 

CO 

SMGOTHEO SPECTRUM WITH INDICATED CONFIDENCE LIMITS 









90 PCT. 




90 PCT. 


95 PCT. 




95 PCT 




COVAR l ANCE 


RAW 


LOWER 




UPPER 


LOWER 




UPPER 


LAG 


SPEC TRUM 


SPECTRUM 


LIMIT 


average 


LIMIT 


limit 


AVERAGE 


LIMIT 


0 


0.9378 


0.3632 


0.3367 


0.3725 


0.4008 


0.3225 


0.3725 


0.4037 


1 


0.6003 


0. 381 7 


0.3256 


0.3602 


0.38 76 


0.3119 


0.3602 


0. 3904 


2 


0.2076 


0.3141 


0.3395 


0. 3755 


0.404 l 


0.3252 


0.3755 


0.4071 


3 


O.OltO 


0.4923 


0.3657 


0.4045 


0.4352 


0.3503 


0.4045 


0.4385 


4 


-0.0592 


0.1191 


0.3115 


0.3446 


0.370B 


0.2984 


0.3446 


0.3735 


5 


-0.04C6 


0.2473 


0.2520 


0.27B6 


0.2999 


0.2414 


0.2788 


0. 3022 


6 


0.0218 


,0. 3011 


0.2537 


0.2807 


0.3020 


0.2431 


0.28C7 


0. 3042 


7 


0.0 79 7 


0. 2732 


0.2560 


0.2B32 


0.3047 


0.2453 


0.2832 


0. 3070 


B 


O.OBCO 


0.2854 


0.2612 


0.28R9 


0.310B 


• 0.2502 


0.2889 


0.3132 


9 


0.0265 


0.3115 


0.2B86 


0.3193 


0. 3435 


0.2765 


0.3193 


0.3461 


10 


-0.0284 


0. 36B6 


0.3366 


0.3723 


0.4006 


0.3224 


0.1723 


0.4036 


1 1 


-0.0670 


0.4406 


0. 3343 


0.3698 


0.3979 


0.3202 


0.1698 


0.4009 


12 


-0.0634 


0.2293 


0.2334 


0.25P2 


0.2778 


0.2236 


0.2582 


0.2799 


13 


-0.0366 


0.1336 


0.1479 


0.1636 


0. 1 761 


0.1417 


0.1636 


0.1774 


14 


-0.0014 


0. 1 580 


0.1 246 


0.1378 


0. 1483 


0.1194 


0.1378 


0.1494 


15 


0.01 76 


0.1016 


0.1045 


0.1156 


0. 1244 


0. IC01 


0.1156 


0.1253 


16 


0.0163 


0. 1013 


0.0926 


0.1025 


0. 1103 


0.0887 


0.1025 


0.1111 


17 


-0.0261 


0.1057 


0 .0635 


0.0923 


0.0994 


0.0800 


0.0923 


0. 1001 


16 


-0.0426 


0.0566 


0.0692 


0.0765 


0.0824 


0.0663 


0.0765 


0.0830 


19 


-0.0340 


0 . OB 7 2 


0.0739 


0.0817 


0.0879 


0.07C8 


0.0817 


0.0886 


20 


-0. 00 7B 


0.0959 


0.0746 


O.0B23 


0.0886 


0.0713 


0.0823 


0.0892 


21 


0.0252 


0.0503 


0.0660 


0.0731 


0.0786 


0.0633 


0.0731 


0.0792 


22 


0.0457 


0.0958 


0.0636 


0.0703 


0.0757 


0.0609 


0.07C3 


0.0762 


23 


0.0485 


0.0394 


0.0544 


0.0602 


0.0648 


0.0521 


0.0602 


0.0653 


24 


0.0337 


0.0662 


0.0529 


0.05B5 


0.0629 


0.0506 


0.0585 


0.0634 


25 


0.016B 


0.0620 


0.0552 


0 • 06 1 0 


0.0657 


0.0529 


0.0610 


0.0662 


26 


0.0096 


0.053B 


0.04BI 


0.0532 


0.0573 


0.0461 


0.0532 


0.0577 


27 


0.0069 


0.0432 


0.0411 


0.0477 


0.05 1 3 


0.0413 


0.04 77 


0.0517 


28 


0.0445 


0.0504 


0.0426 


0.0471 


0.0507 


0.0408 


0.0471 


0.0511 


29 


0.0595 


0.0446 


0.0429 


0.0474 


0.0510 


0.0411 


0.04 74 


0.0514 


30 


0.0164 


0.0502 


0.045R 


0.0506 


0.0545 


0.0438 


0.0506 


0.0549 


31 


-0.0203 


0.0575 


0.0465 


0.0514 


0.0553 


0.0445 


0.0514 


0.0557 


32 


-0.0226 


0.0403 


0.0415 


0.0459 


0.0494 


0.0397 


0.0459 


0.0497 


33 


-0.0139 


0.0453 


0.0409 


0.0452 


0.048 7 


0.0392 


0.0452 


0.0490 


34 


-0.02C4 


0.0500 


0.0425 


0.0470 


0.0506 


0.0407 


0.0470 


0.0510 


35 


-0.0295 


0.0429 


0.0411 


0.0455 


0.0489 


0.0394 


0.0455 


0.0493 


36 


-0.0031 


0.0462 


0.0386 


0.0427 


0.0459 


0.0369 


0.0427 


0.0462 


37 


-0.01 34 


0.0354 


0.0367 


0.0405 


0.0436 


0.0351 


0.0405 


0.0440 


38 


-0.0295 


0.0452 


0.0369 


0.0408 


0.0439 


0.0354 


0.0408 


0.0443 


39 


-0.0343 


0.0375 


0.039B 


0.0441 


0.0474 


0.0382 


0.0441 


0.0478 


40 


-0.0150 


0.0560 


0.0423 


0.0468 


0.0503 


0.0405 


0.0468 


0.0507 
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POWER SPECTRUM ESTIMATION 
t FOR DlSCHf Tf • TIME SERIES 
USING TUKEY-HANNINC WEIGHTING FACTORS 



3MM BEADS AO RCWS MIXED LAG « AO 

ft) 

SMCCTHED POWER SPECTRUM 
WITH 90 PCT . CONFIDENCE LIMITS 
01 234567990 0 



SMOOTHED POWER SPECTRUM 
WITH 95 PCT. CONFIOENCE LIMITS 
2 3 4 5 6 7 8 9 0 
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